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1.0  INTRODUCTION 

The  analysis  of  flow  fields  associated  with  hypersonic  and  other  types  of  reacting  flows 
has  provided  the  impetus  for  the  development  of  computational  methods  applicable  to 
nonequilibrium  flow  fields.  Flows  of  this  nature  are  often  associated  with  hypersonic  aircraft, 
planetary  reentry  vehicles,  arc  heaters,  and  propulsion  systems.  As  Tirsky  (Ref.  1)  has  pointed 
out,  the  flow  fields  associated  with  certain  hypersonic  configurations  (e.g.,  reentry  vehicles) 
experience  both  chemical  and  thermal  nonequilibrium.  Numerical  simulations  of  chemically 
and  thermally  relaxing  flows  require  the  input  of  large  amounts  of  property  data  and  of 
assumed  relaxation  models  (e.g.,  reaction  rates)  which  must  then  be  manipulated  into  a  form 
which  is  accessible  to  a  flow  solver.  While  these  data  are  available  for  most  of  the  materials 
of  concern  at  the  AEDC,  they  are  widely  scattered  and  often  are  presented  in  different 
dimensional  um'ts  and  unfamiliar  forms.  Comparing  the  results  of  codes  which  assume  different 
materia]  properties  and  relaxation  models  is  nearly  impossible.  As  a  result,  the  assembly  of 
diverse  nonequilibrium  data  and  thermochemical  models  into  a  single  package  with  a  common 
format  is  an  essential  part  of  the  computational  fluid  dynamics  (CFD)  effort  at  the  AEDC. 

Since  the  use  of  such  a  common  database  fixes  most  of  the  input  formats  for  the  user 
code  and  essentially  determines  the  methods  for  calculating  transport  and  thermal  properties, 
efflcient  routines  for  performing  these  functions  can  be  prepared  in  advance  and  associated 
with  the  database  in  a  “chemistry  package.”  Furthermore,  since  there  is  virtual  unanimity 
on  the  functional  form  of  the  reaction  rate  —  at  least  for  elementary  reactions  —  this  package 
can  also  contain  a  subroutine  for  evaluating  the  source  terms  in  the  species  conservation  and 
energy  conservation  equations. 

A  number  of  such  “chemistry  packages”  (Refs.  2-4)  is  available.  The  particular  package, 
NEQPAK,  described  in  this  report  is  unique  because  it  allows  for  thermal  nonequilibrium, 
for  muUiple^temperature  reaction  models,  and  for  very  high  temperatures.  NEQPAK  consists 
of  a  database  of  reaction  rates  and  molecular  properties,  and  subroutines  to  evaluate  the 
properties  as  required.  In  addition,  there  are  subroutines  to  evaluate  the  source  terms  and 
their  derivatives  as  well  as  utility  routines  to  convert  units  and  to  compute  some  frequently 
used  functions  of  the  composition.  In  essence,  NEQPAK  provides  the  facilities  needed  to 
evaluate  the  right-hand  sides  of  the  reacting  flow  equations  without  limiting  the  choice  of 
flow  solver  algorithms. 

A  preliminary  version  of  NEQPAK  has  been  incorporated  into  a  one-dimensional  flow 
solver  (Ref.  S)  which  has  been  shown  to  compare  very  well  with  results  obtained  with  other 
codes  (Ref.  6). 
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The  aerophysical  models  embedded  in  NEQPAK  are  discussed  in  Sec.  2.0,  These  models 
have  been  extended  to  include  the  effects  of  ambipolar  diffusion  and  nonequilibrium  thermal 
conductivity  in  Appendixes  A  and  B.  The  extent  and  validity  of  the  database  are  discussed 
in  Sec.  2.7.  An  overall  picture  of  NEQPAK  and  how  it  might  fit  into  a  typical  CFD  code 
may  be  obtained  by  following  the  “scenario”  of  Appendix  I.  Detailed  information  required 
for  inclusion  of  NEQPAK  routines  into  a  flow  solver  is  presented  in  Secs.  3. 1-3. 9. 
Recommendations  for  continuing  the  availability  and  utility  of  NEQPAK  as  a  state-of-the- 
art  tool  are  given  in  Sec.  4.0.  The  Fortran  variables  in  common  are  listed  in  Appendix  C. 
All  NEQPAK  modules  are  listed  in  Appendix  D.  Appendix  E  discusses  the  procedure  for 
providing  the  reaction  model  input.  Sample  reaction  input  sets  are  listed  in  Appendix  F.  These 
reactions  are  not  a  part  of  NEQPAK  per  se,  but  are  typical  of  the  reaction  data  which  the 
user  must  provide  to  any  program  which  uses  the  NEQPAK  package.  Appendix  G  presents 
the  material  database  format.  The  species  for  which  data  are  currently  available  are  listed 
in  Appendix  H. 


2.0  AEROPHYSICAL  MODELS 

NEQPAK  is  designed  to  aid  in  the  analysis  of  flow  of  a  reacting  mixture  of  perfect  gases, 
i.e.,  gases  which  may  be  described  by  the  thermal  state  equation. 


Pi  =  isKT 


(1) 


where  7,  is  the  concentration  of  species  5  in  kg-mole/m^  and  =  8314.34  J/kg-mole/K  is 
the  universal  gas  constant.*  Dalton’s  law  of  partial  pressures 


Na 

P= 

4=1 


(2) 


where  is  the  number  of  species,  is  also  assumed.  Equations  (1)  and  (2)  do  not  adequately 
describe  the  behavior  of  saturated  vapors  or  of  gases  under  very  high  pressures  and 
low  temperatures. 


2.1  REACTING  FLOW  EQUATIONS 


The  equations  which  describe  the  flow  of  a  reacting  compressible  gas  may  be  derived 
from  the  principles  of  kinetic  theory  (Ref.  7).  These  equations  provide  for  the  conservation 
of  species:  ^ 

^  +  Vi(7i£^‘)  =  -Vi(7i«‘)  +  wa,  a  =  1, .. . ,  JV. 


*  SI  units  are  used  throughout  this  report  and  in  NEQPAK  itself. 
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momentum; 

+  Vi(pU*U^)  = 

+  +  +  i  =  1,2,3 

J=1 

and  energy: 


dt 


(4) 


djpE) 

dt 


+  ViipHW)  = 


+Vj[vuH^jV^  +  ViW)  -  IvSijU'VkU*^]  -  g. 


(5) 


where  V,-  =  d/dx*  and  the  summation  convention  is  employed.  The  concentration  of  species 
s  is  denoted  by  and  the  rate  of  change  of  concentration  due  to  chemical  reactions  by  Wj 
in  conformity  with  the  customary  notation.  Thermodynamic  quantities  such  as  the  species 
enthalpy  hf  and  the  total  enthalpy  of  the  mixture  H  usually  have  the  dimensions  J/kg-mole; 

A  * 

mass-specific  values  are  denoted  by  the  symbol  (),  e.g.,  A.  The  mass-averaged  diffusion  velocity 
of  the  s-th  species  is  denoted  by  The  body  force  which  might  represent  gravitational 
forces  or  an  externally  applied  electromagnetic  field,  is  often  neglected.  The  viscosity  of  the 
gas  mixture  is  represented  by  tj,  and  the  thermai  conductivity  of  the  fth  energy  mode  by 
The  number  of  nonequilibrium  energy  modes  is  denoted  by  N/. 


The  conservation  equations  stated  above  are  quite  general;  they  are  limited  by  the 
assumptions  that  the  medium  is  a  continuum  and  that  it  is  possible  to  define  a  translational 
temperature  T  which  is  locally  isotropic.  In  addition,  it  is  assumed  that  the  molecules  possess 
internal  degrees  of  freedom  and  associated  internai  energies  which  may  be  characterized  by 
Boltzmann  temperatures  These  internal  temperatures  may  differ  from  each  other  and 
from  the  translational  temperature  T.  If  thermal  nonequilibrium  is  to  be  simulated,  a  separate 
conservation  equation  (Refs.  8-9) 

+  V.(pe'tf*)  =  Vi(VVr,)  -  Vi(£  -I- 

.=1  (6) 
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with  a  source  term,  uf  =  cIjj  ,  must  be  added  for  each  assumed  nonequilibrium  internal 
energy  mode.  The  limits  m  and  n  will,  in  general,  depend  on  /.  Each  source  term  may  be 
a  sum  of  several  effects.  For  instance,  if  vibration  can  be  described  as  a  Boltzmann  distribution 
with  a  parameter  ^  T,  the  vibrational  energy  equation  for  species  s  would  be 


dt 


+  =  Vi(A*'VT,.,)  - 


(7) 


where  the  source  term  might  be 

•,u  _  c'^-T  1  c^—v  I  ffv—e  1  c"— c 

w, -5,  +S,  +S,  +5,  (8) 

where  5}“^  represents  the  interaction  with  the  translational  energy,  which  is  modeled  in 
NEQPAK  with  the  Landau-Teller  theory,  S""’'  represents  the  the  resonant  interchange  of 
vibrational  energy  with  other  vibrators,  and  is  the  energy  gain  by  electron  Impact.  The 
Final  term  of  Eq.  (B),  5^““^,  represents  the  effect  of  chemical  changes  affecting  the  population 
of  species  s.  Since  dissociation  of  vibrationally  excited  molecules  is  more  likely  and 
recombination  probably  leaves  the  product  molecule  in  an  excited  state,  the  average  amount 
of  vibrational  energy  transferred  by  chemical  interaction  is  not  clear  at  the  present  time. 
NEQPAK,  therefore,  offers  the  user  a  choice  among  the  formulae  advocated  in  Refs.  10-12. 
If  the  vibrational  energy  has  a  non-Boltzmann  distribution,  the  last  term  of  Eq.  (8),  as 
programmed  in  NEQPAK,  should  be  replaced,  but  the  other  three  terms  would  still  apply. 
It  is  frequently  necessary  to  provide  a  conservation  equation  with  source  term. 


St- 


1,  sr 

molecules 


(9) 


for  the  translational  energy  of  the  free  electron,  which  is  affected  by  collision  with  heavy 
particles  by  impact  ionization  S^~^,  and  by  excitation  of  the  vibrational  levels  of  the 
heavy  species,  L^g/ecufes  5}'®.  Both  vibrational  and  electronic  excitation  energy  are  affected 
by  the  emission  and  absorption  of  radiant  energy,  but  these  effects  are  not  fully  implemented 
in  NEQPAK. 


The  foregoing  source  terms  are  described  in  more  detail  in  Secs.  2.5  and  2.6.  The 
implementation  of  these  source  terms  within  NEQPAK  is  described  in  Secs.  3.4  and  3.5. 

The  conservation  equations  must  be  supplemented  by  the  thermal  equation  of  state,  Eq. 
(1),  and  a  caloric  state  equation, 

h,  =  hs(T,  T/)  =  e,  +  1ZT  (10) 

for  each  species.  Since  the  thermal  energy  is  a  nonlinear  function  of  the  temperature  for 
all  species  except  the  electron,  it  is  supplied  by  NEQPAK. 
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The  species  conservation  equations  are  often  written  in  terms  of  species  densities  g, 
=  where is  the  molecular  weight  of  species  s.  Another  common  composition 

variable,  useful  where  the  flow  solver  does  not  require  conservation  form,  is  7^  =  Y/g  where 
g  =  E^igj  is  the  mixture  density.  Mass  fractions  =  Q/q  and  mole  fractions  Xj  =  Y/E^ 
Y^  are  sometimes  employed;  they  have  the  advantage  of  displaying  chemical  processes 
independently  of  the  density  changes  caused  by  the  flow.  Note  that  the  conserved  variable, 
whatever  the  composition  or  energy  units,  is  always  a  density:  number  density,  mass  density, 
or  energy  density. 

The  flow  equations  have  been  introduced  here  to  define  the  notation  and  to  supply  a 
framework  for  the  discussion  of  NEQPAK  capabilities.  NEQPAK  does  not  undertake  to 
solve  these  equations;  it  is  a  set  of  tools  designed  to  save  the  code  developer  time  and  effort 
by  eliminating  the  requirement  for  developing  thermal  and  chemical  routines. 

2.1.1  Boundary  Conditions 

The  conservation  and  state  equations  must,  of  course,  be  supplemented  by  boundary 
conditions  for  each  variable.  Even  if  one  neglects  injection,  ablation,  and  surface  combustion, 
the  ability  of  a  solid  surface  to  absorb  energy  and  to  accumulate  atoms  and  free  radicals 
profoundly  affects  the  chemical  processes  in  its  immediate  vicinity.  One  common  treatment 
of  these  wall  catalytic  effeas  is  to  assume  that  the  gas  at  the  wall  is  in  complete  thermochemical 
equilibrium  at  the  local  wall  temperature  and  pressure.  NEQPAK  has  the  capability,  described 
in  Sec.  3.2.2,  of  generating  a  table  of  equilibrium  compositions  corresponding  to  specified 
pressures  and  temperatures. 

2.2  THERMAL  MODELS 

The  number  and  complexity  of  the  internal  energy  conservation  equations  have  stimulated 
the  formulation  of  numerous  approximate  models.  For  instance,  Lee’s  thermal  model 
(Ref.  8)  assumes  that  the  vibrational  energies  of  all  diatomic  air  species  may  be  categorized 
by  the  same  Boltzmann  temperature  and  that  the  Boltzmann  temperature  Tg,  which 
describes  the  populations  of  the  bound  electronic  states  of  all  of  the  heavy  particles,  is  assumed 
to  be  the  same  as  the  translational  temperature  of  the  free  electrons.  Park  (Refs.  13-14)  makes 
the  further  simplification  that  Tg  =  s  On  the  other  hand,  Candler  and  MacCormack 
(Ref.  15)  have  sought  a  greater  realism  by  assigning  separate  vibrational  temperatures  to  some 
of  the  species.  All  of  these  models  assume  that  the  rotational  degrees  of  freedom  are  in 
equilibrium  with  the  translational  temperature  and  that  the  separate  energy  storage  modes 
are  independent.  Liu  and  Vinokur  (Ref.  16)  argue  that  the  known  dependence  of  the  vibrational 
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potential  energy  upon  both  the  rotational  and  electronic  state  indicates  that  only  a  single 
internal  energy  mode  should  be  considered.  The  following  thermal  models  are  embodied 
in  NEQPAK: 


•  Model  1  corresponds  to  complete  thermal  equilibrium.  All  of  the  internal  energy 
modes  of  all  species  are  in  equilibrium  with  the  translational  temperature  T.  Thus, 
only  one  energy  equation,  Eq.  (S),  is  needed. 


•  Model  2  corresponds  to  Park's  two-temperature  model  (Refs.  13-14)  as  discussed 
in  Ref.  9.  The  translational  and  rotational  modes  are  assumed  to  be  in  equilibrium 
with  the  translational  temperature.  The  electronic  excitation  energy  of  all  species 
and  the  vibrational  energies  of  the  polyatomic  species  are  calculated  with  Ty  as 
is  the  translational  energy  of  the  free  electrons.  Thus, 


where  A®  ^  is  the  molar  heat  of  formation  at  0  K  and  1  atm,  is  the  (constant) 
heat  capacity  of  the  translation-rotation  mode,  and  Cp^  is  the  molar  heat  capacity 
of  the  vibration  and  electronic  excitation  modes.  Note  that  the  heat  of  formation 
is  assigned  to  the  term.  For  species  other  than  the  electron,  the  “pV  work” 
term 72 r  is  assigned  to  h''.  In  the  case  of  the  free  electron,  the  Xtrm'KTg  is 
grouped  with  hg  for  Model  2  and  with  for  Models  3  and  4.  In  addition  to 
Eq.  (5),  a  conservation  equation  for  the  mixture  internal  energy  j 

must  be  supplied.  The  source  term  cL'*'  is  formed  by  summing  the  vibrational 
source  terms,  Eq.  (8),  over  the  molecules  and  adding  the  free  electron  source 
term,  Eq.  (9).  In  the  summation,  the  terms  cancel  out. 


■  Model  3  corresponds  to  the  three-temperature  model  suggested  by  Lee  (Ref.  8) 
in  which  the  distribution  of  vibrational  energy  among  the  available  states  for  each 
molecule  is  described  by  the  single  Boltzmann  temperature  and  the 
populations  of  the  excited  electronic  states  of  all  species  are  described  by  the 
Boltzmann  temperature  7^.  The  kinetic  energy  of  the  free  eleetron  is  assumed 
to  be  in  equilibrium  with  this  electronic  temperature.  Thus, 


A. 


=  h^:  + 


rTv  fTt 

/  c;_gdr+  ci,dr 

yO  .  y® 


(12) 


where  ^  and  Cp  ^  are  the  corresponding  heat  capacities.  Three  energy  equations, 
one  for  the  total  energy,  Eq.  (5),  one  for  the  mixture  vibrational  energy 
and  one  for  the  electron-electronic  energy,  are  needed. 
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•  Model  4  corresponds  to  the  model  employed  by  Candler  and  MacCormack  (Ref. 
15)  in  which  a  separate  vibrational  temperature  ^  is  assigned  to  each  of 
molecules.  The  vibrational  energy  of  the  remaining  molecules  is  assumed  to  be 
in  equilibrium  with  the  translational  temperature.  All  species  share  the  same 
electron  temperature  T^.  There  will  be  one  total  energy  equation,  one  mixture 
electron-electronic  energy  equation,  and  vibrational  energy  equations  with 
source  terms  given  in  Eq.  (8). 

2.3  TRANSPORT  PROPERTIES 


The  species  diffusion  velocity  which  appears  in  Eqs.  (3),  (5),  and  (6)  is  determined  from 


with 


E 


n=l 


XnX« 


("n 


^(7^--^)lvinr, 

Pn  Pa 


(13) 


da  =  Vx,  +  (x*  -  Ti)Vlnji  -  Ya{pT,  -  Y^Pk^kMp  ^14) 

where  is  the  body  force  per  unit  mass,  x^  is  the  mole  fraction  of  species  s  in  the  mixture, 
and  Yg  is  the  mass  fraction.  In  addition,  the  binary  diffusion  coefficients  are  denoted  by 
and  the  thermal  diffusion  coefficients  by  . 


Since  the  solution  of  Eq.  (13)  is  time  consuming,  it  is  customary  to  introduce  the 
approximation  (Ref.  17) 

Xaf^a  ~  —J^adt  (15) 

where 

D  - 


and  dj  is  defined  in  Eq.  (14).  The  quantity  is  called  the  effective  diffusion  coefficient  for 
species  s.  Although  Eq.  (15)  was  derived  under  the  assumption  that  only  a  trace  amount 
of  species  s  is  present  in  the  mixture,  the  diffusion  velocities  calculated  with  Eq.  (15)  agree 
well  with  an  exact  solution  of  Eq.  (13)  for  air  at  7  =  8141 .6  K  as  shown  in  Table  1.  The 
indicated  errors  in  Table  1  are  no  worse  than  the  uncertainties  in  the  binary  diffusion 
coefficients.  This  agreement  may  be  expected  when  the  molecular  masses  of  all  species  are 
of  the  same  order  of  magnitude  (Ref.  18).  If  the  molecular  masses  differ  widely,  use  of  the 
effective  diffusion  coefficients  sometimes  produces  serious  local  errors  in  the  composition 
and  in  mass  conservation  (Ref.  19).  Moss  (Ref.  20)  found  that  the  inclusion  of  the  full 
multicomponent  diffusion  mechanism  made  essentially  no  difference  in  his  high-altitude  viscous 
shock  layer  solutions. 
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Table  1.  Comparison  of  Diffusion  Veiocities  for  High-Temperature  Air 


MolC*¥i| 
Eqs.  (13H14) 
Eqs.  <15~(16) 


1 - 

!  E- 

N 

O 

4.488E-04 

-6.772E-01 

-6.770E-01 

3.225E-OI 
-  2.690E-03 
-2.580E-03 

2.78SE-0I 
-  1.490E-03 
-I.322E-03 

I.405E-02 

7.632E-04 

MlSE-04 


4.488E-04 
-2.176E-03 
-  1 .458E-03 


The  term  ind^,  Eq.  (14),  which  is  proportional  to  V/n  p,  is  customarily  neglected.  The 
thermal  diffusion  term  in  Eq.  (13)  is  usually  insignificant  compared  to  the  species  gradient 
term  and  will  be  neglected  here. 


in  an  ionized  gas,  the  strong  electrostatic  forces  between  positive  and  negative  ions  tend 
to  give  these  particles  a  common  velocity.  This  effect  is  known  as  ambipolar  diffusion 
and  applies  everywhere  except  very  close  to  a  solid  wall.  The  limiting  distance  is  the  Debye 
length,  L 

L  —  (17) 

where  is  the  electron  temperature  in  K  and  is  the  electron  number  density  in  It 
is  shown  in  Appendix  A  that  this  ambipolar  diffusion  introduces  a  further  coupling  between 
the  species  equations,  Eq.  (3),  which  must  now  be  rewritten 


+  V  ■  (7.^.)  =  tii,  V  •  [r(£>,Vx. 


»  V, 


where  is  the  number  of  excess  electrons  in  species  s  and  is  negative  for  a  positive  ion, 

r  =  7^,  and  /i,  is  the  ionic  mobility  of  species  s  and  is  related  to  the  effective  diffusion 

coefficient  by  the  “Einstein  relation”  (Ref.  21) 


fis  =  1 1605.01?, /r. 


(19) 


Note  that  if  only  electrons  and  one  species  of  positive  ions  are  considered,  Eq.  (18) 
reduces  to 

^  +  V .  (y+u)  =  V .  (rz?„yx+) 

where 

+P+Z?, 
fi+  + 

If  it  is  assumed  that  Hg  >  T,  and  Einstein’s  relation,  Eq.  (19),  is  employed,  the 

approximate  result 
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w  2Z)+ 


used  by  Gnoffo,  et  al.  (Ref.  9),  is  recovered. 


Diffusion  and  viscosity  depend  only  on  the  mass  and  motion  of  the  molecules  as  a  whole 
and  are  independent  of  the  internal  state.  Thermal  conductivity,  on  the  other  hand,  is 
determined  by  the  transport  of  energy  —  kinetic  and  internal  —  by  the  molecules. 


Appendix  B  shows  that  the  thermal  conductivity  of  a  pure  species  is  related  to  the  viscosity 
and  molar  heat  capacity  of  the  species  by 


—  [3.7577*  "f 


(22) 


The  first  terra  on  the  right  side  of  Eq.  (22)  is  the  translational  part  X/  of  the  “frozen” 
thermal  conductivity  and  is  unaffected  by  any  thermal  nonequilibrium.  The  second  term, 
on  the  other  hand,  represents  the  contribution  of  the  internal  energy  modes  to  the  thermal 
conductivity  and  is,  therefore,  affected  by  thermal  nonequilibrium  through  Cp^.  As  written, 
Eq.  (22)  is  applicable  to  a  single-temperature  thermal  model.  Note  that  in  this  case, 
represents  the  electronic  excitation  energy  for  both  monatomic  and  polyatomic  species  as 
well  as  the  vibrational  and  rotational  energy  modes.  For  multitemperature  thermal  models, 
a  separate  coefficient 


-  oT>  ^ 


(23) 


is  calculated  for  each  nonequilibrium  mode. 

Note  that  the  transport  properties  of  the  heavy  species  are  evaluated  with  the  translational 
temperature  and  that  77^  and  X^  are  evaluated  with  the  translational  temperature  of  the  free 
electron  Tg-  For  Model  1,  the  user  must  set  Tg  =  T. 

Only  the  viscosity  and  thermal  conductivities  of  the  mixture  as  a  whole  enter  the 
conservation  equations.  The  “exact”  formulae  for  the  viscosity  and  conductivity  of  a  mixture 
in  terms  of  the  species  properties  (Ref.  7)  are  so  unwieldy  and  computationally  expensive 
that  numerous  approximations  have  been  proposed.  The  procedures  adopted  in  NEQPAK 
are  those  of  Armaly  and  Sutton  (Refs.  22-23)  which  are  the  only  approximate  rules  to  explicitly 
account  for  ionization.  The  algorithms  suggested  by  Armaly  and  Sutton  are  too  lengthy  for 
discussion  here. 

Evaluation  of  transport  properties  using  NEQPAK  is  discussed  in  Sec.  3.7. 
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2.4  CHEMICAL  SOURCE  TERMS 

A  typical  reaction  —  the  rth  reaction  —  may  be  written  in  the  form 


^ - 


(24) 


in  which  M,  represents  one  mole  of  species  i  and  the  v^i  and  are  the  stoichiometric 
coefficients.  The  species  on  the  left-hand  side  of  Eq.  (24)  are  called  the  reactants  and  those 
on  the  right-hand  side  are  the  products,  although  this  nomenclature  is  arbitrary  since  the 
reaction  is  also  going  from  right  to  left  as  indicated  by  the  shorthand  symbol  ** .  The  rate 
of  disappearance  of  any  species  involved  in  the  /th  reaction  is 


Ct^  =  n  -I-  vr,,kl  n 

j=i  i=i 


(25) 


A  product  in  reaction  r  appears  at  the  rate 


=  1-;^  ft  7.'''' + -r,.*;  n  li''. 

i=l  i=l 


(26) 


The  net  rate  of  change  of  species  s  due  to  all  reactions  is 

Nr 


=  E(Sr..  -  £...)• 


T-=l 


(27) 


where  Nj.  is  the  number  of  reactions.  Like  most  models  of  the  aerothermochemical  effects 
of  atmospheric  entry  and  combustion,  NEQPAK  assumes  the  modified  Arrhenius  form 
(Ref.  24)  for  chemical  reaction  rates 

hi  =  A,rf'€a:p(-Cr/T,),  (28) 

where  represents  a  generic  temperature.  The  effects  of  thermal  nonequilibrium  on  these 
rates  are  often  modeled  (Refs.  13-14)  by  using  some  blend  of  the  translational  and  internal 
temperatures  for  in  certain  reactions.  Park’s  model  (Ref.  13),  for  example,  uses 
=  T\T^  0  <  j3  <  1 .  The  electron  translational  temperature  Tg  may  be  used  as  for 
electron-impact  ionization  or  dissociative  recombination  reactions.  An  alternative  approach, 
explored  in  detail  in  Sec.  2.4.2,  is  to  assume  that  Eq.  (28)  represents  the  rate  for  thermal 
equilibrium  and  must  be  corrected  by  a  multiplicative  “Vibrational  Coupling  Factor.’’  This 
approach,  proposed  by  Treanor  and  Marrone  (Refs.  25-26),  has  been  extended  recently  by 
Olynick  and  Hassan  (Ref.  1 1)  and  by  Knab,  Friihauf,  and  Jonas  (Ref.  10).  The  ‘virtual  spedes* 
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approach  of  Landrum  and  Candler  (Ref.  12)  results  in  rate  equations  which  can  be  manipulated 
into  the  coupling  factor  form. 

2.4.1  Calculation  of  Reverse  Rates 


(29) 

(30) 


(31) 


In  Eq.  (30),  Pa„„  =  1.01325  x  10^  Pa.  Also,  Cf  is  the  standard -state  Gibbs  free  energy  for 
species  s. 

Many  collections  of  reaction  rate  data  provide  explicit  expressions  of  the  Arrhenius  form 
for  the  reverse  rate  coefficient.  Explicit  specification  of  both  forward  and  reverse  rates 
determines  an  equilibrium  constant,  through  Eq.  (29),  and  a  Gibbs  free  energy,  Eq.  (30), 
which  in  general  will  differ  from  that  calculated  with  the  NEQPAK  curve  fits  and  may  affect 
the  calculated  state  of  the  gas. 

2.4.2  Vibrational  Coupling 

The  vibrational  coupling  method  of  modeling  the  effect  of  vibrational  excitation  on 
chemical  reaction  rates  assumes  that  the  quoted  rate  coefficients  were  measured  for  a  thermal 
equilibrium  situation  and  multiplies  Eq.  (28)  by  a  function  of  translational  and  vibrational 
temperatures.  NEQPAK  extends  this  model  slightly  to  allow  separate  coupling  factors  for 
the  forward  and  reverse  directions. 

Coupling  models  proposed  by  Olynick  and  Hassan  (Ref.  11),  by  Landrum  and  Candler 
(Ref.  12),  and  by  Knab,  Fruhauf,  and  Jonas  (Ref.  10)  are  implemented  in  NEQPAK.  The 
foregoing  coupling  models  are  mathematically  consistent  with  the  use  of  a  common  vibrational 
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temperature  (Model  2  or  Model  3),  but  the  use  of  such  a  crude  thermodynamic  approximation 
largely  negates  any  advantage  which  might  be  obtained  from  the  use  of  these  sophisticated 
coupling  models.  The  Landrum-Candler  and  Knab-Friihauf-Jonas  coupling  models  are 
implemented  in  NEQPAK  only  for  complete  vibrational  nonequilibrium  (Model  4). 


Olynick — Hassan  Coupling  The  Olynick-Hassan  coupling  model  applies  only  to  dissocia¬ 
tion  reactions 

where  molecule  AB  is  a  rotationless  harmonic  oscillator  and  M  is  an  arbitrary  third  body. 
Since  the  analysis  does  not  allow  for  population  inversions,  the  Olynick-Hassan  model  is 
not  valid  for  expanding  flows. 


Let  z  =  and  define 

x  =  {\-TJT)T^J% 

where  ^  is  the  dissociation  temperature  of  AB.  Then,  the  forward  coupling  factor  is 

.p.-  ^r(z,ar) 

^  **  (32) 

where  is  the  incomplete  gamma  function  (Ref.  27).  The  reverse  coupling  factor  is 

=  1-0 

Landrum-Candler  The  Landrum-Candler  model  was  developed  from  an  intensive  theoretical 
investigation  of  the  dissociation  of  Nj.  Its  applicability  to  other  diatomic  species  has  not 
been  investigated,  so  the  generalization  coded  here  may  not  be  warranted.  The  essence  of 
the  model  is  to  replace  the  ordinary  dissociation  by  a  two-step  process 


AB^-M  ^  AB*  +  M  (33) 

AB*  +  M  ^  A  +  B  +  M 


In  Eq.  (33),  AB*  represents  the  AB  molecules  with  vibrational  quantum  numbers  v  >  v^, 
=  {Tfi  AB  ~  where  T^^ab  -  dissociation  temperature.  The  molecules 

AB  have  a  vibrational  temperature  ^  T,  but  the  virtual  species  AB*  is  in  thermal 
equilibrium.  It  follows  that  a  vibrational  coupling  factor  is  needed  only  for  AB.  The  coupling 
factor  for  the  first  of  the  reactions  listed  in  Eq.  (33)  is 


=  exp(/9fla 


T  —  rD.exp[— (Ue 
TyT 


(34) 
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In  Eq.  (34),Z^ff  is  the  partition  function  for  a  harmonic  oscillator  truncated  at  a  quantum 
number  =  [e/(?25^g)]  -  1  below  some  cut-off  energy  e,  namely 


^AB(x,y) 


ImAX  1  j/u 

y-  e-ioAB/v  = 

“  1  -  €-^ABh 


(35) 


where  x  =  ^ab  **  characteristic  temperature  for  oscillator  AB,  and  >'  is  a 

temperature-like  function.  The  parameter  is  supposed  to  correct  for  the  anharmonicity 
of  >15.  Landrum  and  Candler  were  able  to  obtain  impressive  agreement  with  the  full  simulation 
results  for  N2  by  taking  5  =  5.0  and  increasing  the  Millikan-White  (Ref.  28)  vibrational 
relaxation  time  t  by  a  factor  of  25. 


Knab,  Friihauf,  and  Jonas  Knab,  Friihauf,  and  Jonas  (Ref.  10)  have  extended  the  concept 
of  vibrational  coupling  to  reaction  types  other  than  dissociation.  Table  2  shows  for  each 
type  the  coupling  factors  (T.T^)  recommended  by  Knab,  et  al. 


Table  2.  Coupling  Factors  for  Various  Reaction  Types 


Reaction 

Forward  Rate 

Reverse  Rate 

AB  +  M  ^  A  +  B  +  M 

K 

AB  +  C  ^  AC  +  B 

A  +  B  r^AB*  -H  E~ 

In  Table  2,  Ar-f  and  Arj:  are  the  ordinary  equilibrium  reaction  rates  in  the  Arrhenius  form. 
The  ^  functions  are  defined  as  follows: 

^  _  2AB(Td^AB,T)  Zab{0‘t-,9)  -(•  ZAB{Td,AB,To,r)  "  (^r , Tp.r ) 

^  2>tB(Tj,XB»T„)  e-'^^l'^ZAB{f^T,-^r)-\-ZAB{Td,AByT*)-ZAB{*iT^T*) 

^  ^  ZAB{Ti,AB,T)  .. 

^  ZAB{Ti^AB ,  T„^Ab) 

ZAB(<^Td,ABf 9Ab)  +  ZAB{Td,AB^T(i,.r)  -  ZAB{f^’I'd,AB,To^T) 
ZAB{aTd,AB^  “<^r)  +  ■2yifl(rj,^fl,T*)  -  ZAB{oZ'd^AB,'J^)  (37) 
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^  ^0,r) 

^  2,(Trf.„r„,,)3,(rj,„r*)  (38) 

where  s  represents  the  molecule  or  ion  on  the  right  side  of  the  reaction.  The  pseudo- 
temperatures  Tq  ^,  and  T*  used  in  Eqs.  (36)-(38)  are 

J-  =  1  1  1 

9r  r,,,  T  <f>r  (39) 

_1_  _  J _ 1_ 

Tojr  ~  T„,a  4>t  (40) 


J_  _  j__J_ 

T*  ~  T  (i>r  (41) 

The  input  parameter  measures  the  probability  that  a  molecule  in  vibrational  state  / 
will  dissociate.  If  =  *.  dissociation  is  equally  probable  from  all  vibrational  levels.  Knab, 
et  al.  also  introduce  an  effective  activation  temperature  a,  =  aC^  for  reaction  r.  [f  a  =  1, 
the  formulation  of  Eq.  (37)  reduces  to  that  of  Treanor  and  Marrone  (Ref.  25). 

For  associative  ionization  reaaions,  the  vibrational  coupling  factor  for  the  reverse  direaion 
is  just  the  ^  j  function  defined  in  Eq.  (38)  above  with  the  translational  temperature  T  replaced 
with  the  electron  temperature  T^. 

2.5  VIBRATIONAL  ENERGY  SOURCE  TERMS 

The  species  vibrational  energy  is  affected  by  inelastic  collisions  with  heavy  molecules  and 
with  electrons.  During  a  heavy  particle  collision  there  may  be  an  exchange  of  vibrational 
energy.  The  species  vibrational  energy  is  also  affected  by  an  increase  or  decrease  in  the  number 
of  vibrationally  excited  molecules  through  chemical  reaction.  Interaction  with  the 
electromagnetic  field  is  not  modeled  at  the  present  time. 

2.5.1  V-T  Interaction 


Landau  and  Teller  (Ref.  29)  showed  that  if  the  vibrational  energy  levels  are  equally  spaced 
(harmonic  oscillator)  and  only  single-level  transitions  are  allowed,  then 


(42) 
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where  the  indicates  the  equilibrium  value  of  the  vibrational  energy,  and  is  the  so-called 
relaxation  time.  Strictly  speaking,  Eq.  (42)  is  valid  only  for  a  harmonic  oscillator,  but  it  is 
applied  to  anharraonic  oscillators  by  using  the  polynomial  curve  fits  defined  in  Sec.  2.7  to 
calculate  the  equilibrium  vibrational  energy 


The  relaxation  time  r^j  of  species  s  in  a  bath  of  species  J  is  calculated  using  the  expression 


pTsj  =  exp(r;j/T<0 


(43) 


where  p  is  the  pressure.  The  relaxation  time  employed  in  Eq.  (42)  is 


where 


JV. 


is  Park’s  high-temperature  correction  (Ref.  14  ),  and 


C  =  =  2.898677  x  lO"*^ 


(44) 


(45) 


in  SI  units.  The  value  a,  =  10  ®  for  the  effective  cross  section  is  that  suggested  by 
Gnoffo,  et  al.  (Ref.  9). 


The  necessary  coefficients  for  a  particular  interaction  pair  may  be  generated  by  making 
use  of  the  Millikan  and  White  correlations  (Ref.  28): 


pnj  =  -  18.42] 

and 

Aij  =  o.ooiie^iV'dV® 


(46) 

(47) 


where  p  is  the  pressure  in  atmospheres  and  fi  is  the  reduced  mass  of  the  mixture.  In  Eq.  (47), 
6f  =  hvj/x,  where  p,  is  the  characteristic  vibrational  frequency,  h  is  Planck’s  constant,  and 
X  is  the  Boltzmann  constant.  Also,  is  the  reduced  mass  of  the  interacting  molecules  i 
and  j. 
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2.5.2  V-E  Interaction 


An  encounter  between  a  free  electron  and  a  molecule  may  result  in  the  transfer  of  kinetic 
energy  from  the  electron  to  a  vibrational  mode  of  the  heavy  particle.  The  probability  for 
such  an  interaction  is  an  extremely  irregular  function  of  the  electron  kinetic  energy. 
Nevertheless,  Lee  (Refs,  8  and  30)  has  been  able  to  model  the  process  with  an  equation  of 
the  Landau-Teller  type,  [see  Eq.  (42)],  namely 


er 


(48) 


where  ef*  is  the  vibrational  energy  evaluated  at  the  electron  temperature  and  the  relaxation 
time  is 


_ 2 _ 

ne(l  - 


(49) 


where  is  the  electron  number  density,  0^  is  the  characteristic  vibrational  temperature  of 
species  s,  and  kg  is  the  rate  constant  for  a  transition  from  the  vibrational  ground  state  to 
vibrational  level  v. 


Since  the  vibrational  energy  gained  by  the  molecule  is  lost  by  the  electron,  V-E  interactions 
need  not  be  calculated  for  Model  2,  in  which  vibrational  and  electron  energy  are  pooled. 

2.5.3  V-V  laleraction 


It  is  possible  for  two  molecules  to  exchange  vibrational  energy  during  a  collision,  but 
this  process  is  unlikely  unless  the  molecules  happen  to  have  vibrational  levels  v,  and  Vj  with 
nearly  the  same  energy.  This  vibrational  exchange  may  be  described  as  a  reaction  such  as 

^2(1)  -h  /r20(ooo)  =  ^2(0)  +  ^20(001) 

which  describes  the  interchange  of  vibrational  energy  between  //j  and  one  of  the  modes  of 
//jO.  Rate  data  for  such  V-V  reactions  are  available  (Ref.  31),  but  considering  level-to-level 
transitions  of  the  kind  indicated  is  computationally  expensive.  Therefore,  NEQPAK  has 
adopted  the  model  of  Candler  and  MacCormack  (Ref.  15)  in  the  form  proposed  by  Knab, 
et  al.  (Ref.  10) 


Sr'  = 

jjii  V 


(50) 
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where 


gfl./T  _  1 


and  Vjj  -  0.01  are  the  exchange  probabilities  averaged  over  all  vibrational  energy  levels. 


2.5.4  V-C  Interactions 


It  i$  clear  that  the  dissociation  of  a  vibrating  molecule  removes  its  vibrational  energy  from 
the  species  pool  and  that  the  formation  of  a  new  molecule  must  add  some  vibrational  energy 
to  the  pool.  The  simplest  assumption  is  that  the  energy  gained  or  lost  during  the  reaction 
is  just  the  average  vibrational  energy  of  the  species,  namely 


(51) 


This  “nonpreferential  dissociation”  assumption  is  appropriate  for  Models  2  and  3  and  also 
for  the  Olynick-Hassan  coupling  (Ref.  11)  discussed  in  Sec.  2.4.2. 


Landrum — Candler  CouplingThe  Landrum-Candler  theory  (Sec.  2.4.2)  assumes  that  the  lower 
vibrational  energy  levels  v  <  relax  according  to  the  Landau-Teller  theory  while  the  upper 
vibrational  levels  (the  virtual  species  AB*)  are  in  thermal  equilibrium.  The  critical  level  is 


Vc 


^AB  -  T 
^AB 


+  1 


where  is  the  dissociation  temperature  and  6y^g  is  the  characteristic  vibrational  temperature 
for  species  s.  The  level  spacing  of  the  virtual  species  is  ^AB*  ~  ^i^AB  where  —  0.4  is  an 
input  parameter.  Thus,  the  rate  of  change  of  vibrational  energy  is 


where 


—  -  ^AB - + 


dt 


Tab 


(52) 


-  £„) 

r  X  y 


(53) 


In  Eq.  (53),  x  =  exp(v^^g/7)  —  1  and  ^  =  1  -  exp(-&^g/T).  The  reaction  r  is 

AB  +  M^  AB‘  +  M 
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The  function  F  is  defined  as 

F  - 

Z{a,T) 

where 

Z{a,T)  =  «<p(-(»,  - 


and  =  {v„  +  1  - 


KFJ  Coupling  Unlike  the  other  coupling  theories  described  above,  in  the  KFJ  theory  (Ref. 
10)  the  vibrational  coupling  factor  depends  on  both  the  translational  and  the  species  vibrational 
temperatures.  Thus,  in  terms  of  the  gain  and  loss  rates  introduced  in  Sec.  2.4, 


5;- =  /(T,T„,.)(S,,, -£,,,) 


(54) 


where Ty  y)  is  the  quantity  defined  by  Eqs.  (48)  and  (49)  of  Ref.  10.  The  input 

parameter  a  is  defined  in  Sec.  2.4.2  as  are  the  pseudo-temperatures  T  =  g^,  U  ^  and 
Tq  ^.  The  source  term  of  Eq.  (54)  applies  only  to  the  reactions  listed  in  Table  2.  The  total 
source  term  for  is  formed  by  summing  over  all  reactions. 


2,6  ELECTRON-ELECTRONIC  ENERGY  SOURCE  TERMS 


The  energy  of  the  free  electron  is  affected  by  elastic  and  inelastic  collisions  with  heavy 
particles,  by  electron  production,  by  interaction  with  applied  and  induced  electric  fields,  and 
by  radiation  processes.  The  energy  of  electronic  excitation  may  be  lost  to  radiation.  Of  these 
effects,  only  those  dealing  with  collisions  and  electron  production  are  included  among  the 
NEQPAK  source  terms  of  Eq.  (6)]. 


Elastic  collisions  between  electrons  and  heavy  particles  cause  an  energy  exchange  described 
by  term  6  in  Eq,  (4)  of  Ref.  9,  namely 


sr^  =  S1lA4,7,(r  -  T.)  .fr 

(55) 


where  is  the  collision  frequency  defined  by 


/eg 


7e 


(56) 
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for  collisions  with  charged  particles  and  ^ 

/«  =  (57) 

for  collisions  with  neutral  particles.  The  constants  ^2*  ^3  defined  as 


6  =  8\/ir/me 


J\fe* 

3(2k)1-6 


(58) 


6  =  K^/(irA/e®) 


(59) 


^  =  MyJ^Kfivme)  (60) 

where  is  the  mass  of  the  electron,  x  is  Boltzmann’s  constant,  and  e  is  the  elementary 
charge.  The  effective  cross  section  a  is  computed  with  quadratic  fits  taken  from  Ref.  9. 

Electron  impact  ionization,  i.e.,  a  process  such  as 

AB^-E-  — ►  i4J9+  +  £-+!;■ 


leads  to  an  electron  energy  loss 


Sr  S. 


r=l ss2 


(61) 


In  Eq.  (61),  is  the  forward  rate  coefficient  defined  in  Eq.  (28)  andX^  is  the  ionization 
energy  per  kg-mole  of  species  s.  The  product  of  stoichiometric  factors  will  vanish  except 
for  impact  ionization  reactions. 


2.7  THE  MATERIAL  PROPERTIES  DATABASE 


Post-shock  temperatures  in  the  flow  field  of  a  reentering  planetary  vehicle  can  approach 
50, (KK)  K.  While.direct  calorimetric  measurements  at  such  temperatures  are  clearly  impossible, 
it  is  possible  to  use  spectroscopic  data  and  the  methods  of  statistical  mechanics  to  infer  thermal 
properties  at  reentry  conditions.  The  most  extensive  tables  of  thermodynamic  properties  at 
very  high  temperatures  are  those  of  Browne  (Refs.  32-35),  who  used  low-temperature 
measurements  to  evaluate  the  parameters  of  the  Morse  potential  from  which  he  calculated 
the  second  virial  coefficient  and  the  partition  function  (Ref.  36).  This  approach  has  been 
criticized  by  Jaffe  (Ref.  37)  on  the  grounds  that  it  does  not  account  for  the  effect  of  excited 
electronic  states  on  rotational-vibrational  potential  energy  and  that  it  neglects  the  effect  of 
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finite  dissociation  energy  on  the  partition  function.  The  most  accurate  calculations  (limited 
apparently  to  N^)  are  those  of  Ref.  38,  which  involved  thousands  of  numerical  solutions  of 
the  Schrbdinger  wave  equation. 


A  computational  effort  of  the  scale  outlined  above  is  clearly  impracticable  for  all  of  the 
polyatomic  and  ionized  species  in  the  NEQPAK  database.  Moreover,  the  necessary  spectral 
data  are  not  available  for  many  of  the  species  of  interest.  The  JANAF  tables  (Ref.  39)  contain 
tabulated  thermodynamic  data  over  the  range  0  s  T  <  6,000  K  for  over  600  species. 
Polynomial  fits  to  the  JANAF  data  are  presented  by  Svehia,  et  al.  (Ref.  40).  Esch,  et  al. 
(Ref.  41)  have  extended  these  fits  to  15,000  K  for  selected  species  of  special  interest  to  the 
aerospace  community.  Gupta,  et  al.  (Ref.  42)  have  obtained  polynomial  fits  to  Browne’s 
data  up  to  35,000  K  for  II  air  species.  Because  of  the  limited  range  and  inconsistent  quality 
of  the  published  data,  many  of  the  thermodynamic  property  tables  in  NEQPAK  have  been 
recomputed  using  the  methods  of  Ref.  43  with  molecular  constants  taken  from  Refs.  39  and 
44.  The  calculations  were  extended  to  30,000  K  in  order  to  ensure  that  the  polynomial 
approximations  yielded  reasonable  values  at  typical  post-shock  temperatures. 

In  order  to  minimize  memory  requirements,  the  thermal  data  are  stored  in  the  form  of 
multisegment  polynomials 


n 


=  “i,«  +  a2,j,r  -I-  a3,,T*  -f 


(62) 


-  tt,  4.  4-  4.  -L  'T*  A. 

nr  ~  ^  2  ^  ^  3  4  ^  +  T 


(63) 


0i,,(l-lnr)- 


,  06, B, 

20  r  ^ 


“7,* 


(64) 


There  are  up  to  six  segments  for  each  species.  The  junction  temperatures,  which  differ  from 
species  to  species,  mark  the  upper  limit  of  applicability  of  a  given  set  of  coefficients. 


Since  the  convergence  of  iterative  processes  involving  the  thermal  properties  may  be 
hindered  by  discontinuities  in  the  data,  great  care  was  taken  to  maintain  continuity  at  the 
junction  points.  A  special  least-squares  routine  which  enforced  continuity  of  value  and  slope 
at  the  junctions  was  used  to  calculate  the  coefficients  Oj, ...,  a,.  Smoothness  at  the  junction 
points  was  further  improved  by  making  slight  adjustments  in  Og  and  a-,.  The  heat  of 
formation  at  absolute  zero  clearly  dictates  the  value  of  Og  for  the  first  segment. 

Figure  2  is  a  plot  of  the  dimensionless  specific  heat  of  NO  as  calculated  by  Jaffe  (Ref. 
37),  by  Browne  (Ref.  32),  and  with  the  polynomial  fits  used  in  NEQPAK.  The  quality  of 
the  NEQPAK  fits  to  the  enthalpy  may  be  assessed  from  Table  3. 
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Figure  1.  Heat  capacity  of  NO. 


Table  3.  Fit  Quality  for  NO 


Junction 

^upper 

Temperature 

^(ower 

300.00 

1.1508*10-5 

1000.0 

-2.8633*10-* 

5000.0 

2.0833*10-’ 

15000.0 

4.1667*10-® 

25000.0 

5.8333*10-® 

Average  error  0.195% 

Maximum  error  0.581% 

The  database  also  contains  the  coefficients  for  polynomial  Hts  to  the  natural  logarithms 
of  the  binary  diffusion  coefficients  and  viscosities.  These  fits  are  from  various  sources.  The 
fit  coefficients  for  the  eleven  air  species  O2,  N2,  O,  NO,  NO* ,  O2  ,  N2  ,  O* ,  N* ,  and  E~ 
are  taken  directly  from  Ref.  42.  Although  Ref.  42  recommends  the  use  of  the  Sutherland 
formula  for  air  viscosities  at  T  <  1,000  K,  numerical  experiments  at  AEDC  show  that  the 
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tabulated  fits  work  well  at  room  temperature  or  above.  The  viscosity  of  neutral  particles 
other  than  the  air  species  and  the  binary  diffusion  coefficients  involving  these  species  were 
calculated  with  the  program  of  Ref.  19.  Interactions  between  charged  and  neutral  particles 
were  evaluated  with  the  program  of  Ref.  45.  Charged-charged  interactions  were  calculated 
with  a  special  code  based  on  the  formulae  of  Devoto  (Ref.  46). 

The  transport  properties  of  the  individual  species  are  calculated  with  the  following  curve-fit 
expressions  (Ref.  42): 

7*  =  0.1exp(fc5^  -I-  b4^T-\-  1>3^7*  -I-  62,»‘7®  +  bi^7*) 

©ij  =  12:^5  exp(i.,ij  + 

where  T  =  In  T  and  S  is  the  electrostatic  shielding  ratio,  namely 

1  for  neutral-particle  interactions 

i  Ai<ri  either  particle  is  charged 

EnA(r/pc'  ) 

The  shielding  function  A(x)  is  defined  as 

A(a:)  =  (2.09  x  +  1.52  x 

and  Pg  =  1.01325  x  10®  is  the  electron  partial  pressure  expressed  in  atmospheres. 

3.0  USING  NEQPAK 


(65) 

(66) 


3.1  STRUCTURE 

The  data  flow  within  NEQPAK  is  shown  in  Fig.  2.  The  modules  NPCHEMIN  and 
NPREPEQ  shown  in  the  bottom  tier  of  Fig.  2  are  “setup”  routines  called  by  the  flow  solver 
or  a  preprocessor  to  interpret  the  chemical  model  and  to  prepare  the  necessary  constants 
and  tables.  These  setup  routines  need  be  called  only  once  while  modules  in  the  upper  tier 
are  tools  which  are  called  repeatedly  by  the  flow  solver.  Since  it  is  impossible  to  name  all 
of  the  “tool”  routines  in  the  limited  space  of  the  diagram,  only  the  broad  categories  are 
shown.  The  use  of  these  tools  is  explained  fully  in  the  following  sections. 

NPCHEMIN  (Sec.  3.2.1)  reads  the  user-supplied  reaction  mechanism  from  the  file  assigned 
to  unit  number  INPR,  extracts  the  necessary  information  from  the  material  properties  database 
on  unit  numbers  INPTH  and  INREL  and  puts  the  necessary  constants  and  tables  into  the 
common  blocks  /NAMES/  and  /NEQCOM/.  The  thermal  model  and  the  names  and 
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concentrations  of  any  inert  species  are  calling  arguments  of  NPCHEMIN.  In  addition,  if 
Model  4  is  selected,  the  names  of  the  nonequilibrium  vibrators  must  be  specified  by 
the  user. 

The  two  common  blocks  /NAMES/  and  /NEQCOM/  contain  only  the  thermochemical 
and  transport  properties  for  the  current  problem.  The  variables  in  common  block  /NEQCOM/ 
are  defined  in  Appendix  C. 

Subroutine  NPREPEQ  (Sec.  3.2.2)  uses  the  thermochemical  properties  in  the  common 
blocks  to  generate  a  table  of  equilibrium  compositions  which  is  stored  on  unit  lEQ. 
Interpolation  within  this  table  is  performed  by  Subroutine  NPEQCOMP  (Sec.  3.9).  A  third 
common  block  /EQCOM/  contains  information  relative  to  the  number  and  spacing  of  the 
pressures  and  temperatures  for  which  equilibrium  composition  data  are  tabulated. 

All  data  stored  by  the  flow  solver  at  points  of  the  computational  grid,  such  as  local  values 
for  temperature,  energy,  and  composition,  are  communicated  through  the  subroutine  argument 
lists.  In  order  to  save  memory,  these  data  are  passed  to  the  various  NEQPAK  subroutines 
in  arrays  of  vector  length  NP  rather  than  passing  the  data  for  the  entire  grid.  In  a  flow  solver, 
NP  would  ordinarily  be  one  of  the  grid  dimensions.  The  index  referring  to  the  grid  dimension 
is  always  the  last  subscript  of  the  array.  Array  dimensions  are  either  passed  through  the 
subroutine  arguments  or  are  set  by  the  parameter  statement  in  file  NPDIMEN  which  is  included 
in  every  routine.  Table  4  lists  the  current  default  values  for  the  array  dimensions.  Since  some 
flow  solvers  use  a  leap-frog  procedure,  the  arguments  of  all  subroutines  dealing  with  grid- 
dependent  arrays  contain  the  sequence  “LOW,NP,INC  “  in  which  LOW  is  the  index  of  the 
initial  point,  NP  is  the  index  of  the  final  point,  and  INC  is  the  increment  for  “DO-loops.” 

Computation  has  been  vectorized  as  far  as  possible.  This  vectorization  will  be  most  efficient 
when  the  number  of  grid  points  —  always  the  innermost  loop  —  exceeds  the  number  of  species 
or  reactions. 

SI  units  are  used  exclusively  in  NEQPAK.  The  composition  variable  is  the  concentration 
in  kg-moles/m^.  Energies  are  expressed  in  J/kg-mole  or  in  J/m^.  Reaction  rates  may  be 
entered  in  cm^/gfn-mole/s  or  in  cn^/ particle/ s,  but  these  will  be  converted  to  SI  units  before 
being  used  for  calculation. 

Apart  from  error  messages,  only  the  “set-up”  modules,  NPCHEMIN  and  NPREPEQ, 
write  to  standard  output.  NPREPEQ  also  writes  to  unit  lEQ.  Detailed  descriptions  of  the 
various  NEQPAK  routines  are  given  in  the  following  subsections.  Calling  sequences  are 
included  for  these  subroutines  in  the  form  of  tables.  The  tables  give  a  listing  of  the  variables 
in  the  calling  sequences  of  the  subroutines  along  with  the  internal  dimensions  of  any  arrays 
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used  as  input  or  output  by  the  subroutines.  With  one  exception,  these  subroutines  must  be 
called  directly  by  the  flow  solver.  .Although  Subroutine  NPCOUPLE  is  called  only  from 
Subroutine  NTSOURCE,  it  is  considered  here  in  order  to  discuss  the  various  vibrational 
coupling  models  (Sec.  2.4.2)  from  which  the  user  may  choose. 


Material  Properties  Database 
Figure  2.  Schematic  of  NEQPAK  data  flovk. 


Table  4.  Current  Values  of  Some  Important  Dimensions 


MXBR 

Maximum  number  of  thermofit  ranges 

6 

MXNE 

Maximum  number  of  elements 

8 

HXNP 

Maximum  number  of  grid  points  in  a  vector 

101 

MXNR 

Maximum  number  of  reactions 

50 

MXNS 

Maximum  number  of  species 

15 

MXNV 

Ma.ximum  number  of  oscillators 

7 

MXEQT 

Maximum  number  of  equilibrium  temperatures 

26 

HXEQP 

Maximum  number  of  equilibrium  pressures 

26 

LENGTH 

Maximum  record  length  for  unit  lEQ 

8*(MXXS  +  4) 

HXNT 

Maximum  number  of  temperature  types 

9 

HXRAD 

Maximum  number  of  radiating  processes 

0 

HXNI 

Maximum  number  of  internal  energy  modes 

3 
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3.2  SETTING  UP 

In  the  flow  solver,  the  user  must  define  and  open  the  file  which  contains  the  reaction 
mechanism  and  assign  it  to  unit  number  INPR,  while  the  files  which  contain  the  material 
properties  database  must  be  opened  and  assigned  to  unit  numbers  INPTH  and  INREL.  If 
some  region  or  regions  of  the  flow  field  are  assumed  to  be  in  thermochemical  equilibrium, 
the  user  must  also  open  unit  lEQ  for  direct  access  specifying  a  record  length  of  8  *  (MXNS 
+  4)  to  provide  access  to  the  equilibrium  composition  tables  described  below. 

The  user  should  consider  editing  NPDIMEN  to  adjust  the  array  dimensions  for  a  particular 
problem  when  compiling  NEQPAK.  Default  values  for  all  array  dimensions  are  listed  in 
Table  4. 

3.2.1  Input  Processing 

Subroutine  NPCHEMIN  reads  the  reaction  file,  identifies  the  species  involved  in  the 
problem,  and  arranges  them  according  to  category.  The  free  electron  E~ ,  if  present,  is  listed 
first  and  is  followed  sequentially  by  the  other  monatomic  species,  virtual  species  (identified 
by  an  asterisk),  the  NVIB  nonequilibrium  vibrators  listed  in  the  VIB  array,  neutral  molecules, 
and  finally  ions  (positively  or  negatively  charged).  Some  of  these  categories  may  not  exist 
or  there  may  be  some  overlap.  NPCHEMIN  first  identifies  the  species  present  in  the  current 
problem,  and  then  extracts  the  necessary  thermophysical  data  from  the  database,  sets  up 
the  stoichiometric  matrices  relating  the  species  to  the  reactions  in  which  they  appear,  calculates 
certain  species-dependent  constants,  and  sets  all  the  necessary  flags.  These  data  are  loaded 
into  the  common  block  /NEQCOM/.  Certain  items  such  as  species  names  which  might 
be  used  directly  by  the  flow  code  for  reporting  purposes  are  also  returned  through  the 
argument  list. 

Communication  with  NPCHEMIN  is  through  the  calling  sequence  shown  in  Table  5.  Note 
that  VIB,  NVIB,  SPECFS,  and  FSCOMP  are  used  for  both  input  and  output.  Since  some 
of  the  reactants  may  not  have  been  input  to  NPCHEMIN  through  the  SPECFS  array,  the 
number  and  order  of  the  species  may  differ  on  output  from  their  input  values.  The 
compositions  in  FSCOMP  will  be  arranged  to  agree  with  the  species  order  in  SPECFS.  The 
thermal  model  is  selected  by  specifying  a  value  for  the  variable  MODE.  MODE  is  set  equal 
to  1,2,  3,  or  4,  depending  on  which  of  the  four  thermal  models  described  in  Sec.  2.2  the 
user  chooses. 

Input  of  SPECFS  and  FSCOMP  values  is  required  only  for  inert  species  such  as  argon, 
which  do  not  appear  in  the  reactions  being  considered.  Species  names  are  character  constants 
and  must  be  in  standard  chemical  notation  with  no  internal  blanks.  It  is  best  if  the  units 
of  FSCOMP  have  “moles”  or  “particles”  in  the  numerator  such  as  kg-moles/m^. 
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Table  5.  Calling  Sequence  for  Subroutine  NPCHEMIN 


SUBROUTIIE  HPCHBHIS (HALT , IIPR, IHPTH. ZHREL ^NODE ,HSFS .ISPEC , HVZB , 
ITHP . SPECFS , VIB .FSCOKP.UT) 

IIPUT 


HALT 

If  TRUE,  stops  run  if  any  Bissing  data 

IHPB 

Rsaction  input  unit 

IBPTH 

Property  data  input  unit 

IHREL 

Relaxation  data  input  unit 

MODE 

ThexnoeheDical  model 

HSFS 

lumber  of  species  in  fres-stream 

SPECFS(NXHS) 

Fames  of  frse-stream  species 

FSCOMP (NXHS) 

7. 

Ftee-stream  composition 

¥VIB 

lumber  of  nonequilibrium  vibrators 

VlB(MXirV) 

fames  of  nonequilibrlnm  vibrators 

OUTPUT 

5PECFS 

Species  considered  in  this  run 

FSCOMP 

Free  stream  composition  of  above  species 

HSPEC 

Length  of  SPECFS  and  FSCOMP  arrays 

VIB 

fonequilibriuio  vibrators  considered  in  t] 

HVIB 

Length  of  VIB  array 

HTNP 

Hnmber  of  tamperatura  types  considered 

VT(MXNS) 

M, 

Species  molecnlaz  ueights 

Anj  unit  a 


kg/kg-mola 


The  VIB  array  need  only  be  filled  if  MODE  =  4  is  selected.  In  this  case  the  user  must 
set  NVIB  equal  to  the  number  of  vibrators  specified  to  be  in  thermal  nonequilibrium  and 
the  VIB  array  must  contain  the  names  of  the  Ny  vibrators  specified.  On  output  from 
NPCHEMIN,  the  user  must  examine  the  FSCOMP  array  to  determine  the  order  in  which 
NEQPAK  has  indexed  the  vibrators.  If  MODE  =  1 , 2,  or  3  is  chosen,  then  all  the  user  need 
do  is  set  NVIB  =  0  on  input  to  NPCHEMIN. 


Subroutine  NPCHEMIN  will  unconditionally  abort  the  run  if  it  cannot  Find  thermodynamic 
data  for  a  species.  However,  since  NEQPAK  will  run  fairly  well  if  transport  data  are  missing 
for  some  species,  the  logical  variable  HALT  allows  the  user  to  override  program  termination. 
The  run  will  be  terminated  if  any  data  are  missing  and  HALT  =  .TRUE.,  while  the  program 
will  continue  to  run  if  HALT  =  .FALSE. 


In  order  to  avoid  the  necessity  of  testing  at  every  grid  point  at  each  entry  into  some 
subroutines  for  the  very  special  case  that  7,  =  0,  NPCHEMIN  will  add  10“  to  all  elements 
of  the  output  FSCOMP  array.  This  procedure  will  lead  to  an  initial  charge  imbalance  if 
electrons  are  involved  in  the  chemical  model  and  FSCOMP  is  expressed  in  mass  units  such 
as  partial  densities  or  mass  fractions. 
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3.2.2  Equilibrium  Composition  Tabie 

In  reacting  flow  problems,  it  is  necessary  to  augment  the  usual  boundary  conditions  on 
velocity  and  total  energy  with  wall  boundary  conditions  for  each  species.  One  frequent 
assumption  is  that  of  a  “fully  catalytic”  wall,  i.e.,  that  the  gas  is  in  chemical  equilibrium 
at  the  wall  temperature  and  pressure.  This  equilibrium  composition  may  be  calculated  by 
minimizing  the  Gibbs  free  energy  of  the  mixture  subject  to  the  eonstraints  of  a  Hxed  total 
mass  and  fixed  atomic  mass  fractions. 

Since  the  minimization  procedure  is  an  iterative  process  (typically  20  to  30  cycles  per  case) 
involving  the  evaluation  of  the  free  energy  for  each  of  the  n  species  and  the  solution  of  an 
n  by  (n  + 1)  matrix  equation  at  each  cycle,  the  direct  minimization  process  at  a  boundary 
point  is  probably  equivalent  in  CPU  time  to  the  chemical  calculations  at  10  Held  points.  The 
obvious  way  to  overcome  this  delay  is  to  interpolate  in  a  table  which  has  been  specialized 
to  a  small  range  of  temperatures  and  pressures  typical  of  the  problem  in  hand.  For  instance, 
if  wall  temperature  is  specified,  then  pressure  would  be  the  only  independent  variable. 

Subroutine  NPREPEQ  (Table  6)  has  been  set  up  to  create  such  a  table,  although  it  could 
be  used  one  pressure-temperature  point  at  a  time.  NPREPEQ  is  a  driver  program  which 
performs  the  necessary  initialization,  such  as  the  calculation  of  the  atomic  mass  fractions, 
sets  up  the  pressure-temperature  matrix,  calls  Subroutine  NPEQLBRM  to  perform  the  free- 
energy  minimization  for  each  pressure-temperature  point,  and  returns  as  output  the  calculated 
density,  equilibrium  sound  speed,  and  equilibrium  composition.  Only  the  substances  in  the 
species  list  created  by  NPCHEMIN  are  considered  in  the  minimization  procedure.  Species 
free  energies  are  calculated  by  Subroutine  NPGIBBS  and  therefore  are  consistent  with 
NEQPAK  thermodynamics. 

A  convenient  way  of  setting  up  an  aerodynamic  calculation  would  call  NPCHEMIN  and 
NPREPEQ  at  an  early  stage  and  use  the  species  names  output  through  the  NPCHEMIN 
argument  “SPECFS”  to  initialize  the  composition.  A  call  to  NPTHERM  (Sec.  3.6)  could 
then  give  initial  values  for  the  energy  variables. 

3.2.3  Multiple  Domains 

Many  aerodynamic  problems  involve  regions  of  essentially  different  chemistry  such  as 
the  core  of  a  rocket  plume  or  the  boundary  layer  on  an  ablating  or  diffusion-cooled  surface. 
Here,  the  user  might  choose  to  merge  the  applicable  reaction  sets  and  work  with  a  large  species 
list.  Such  a  procedure  would  be  wasteful  of  both  memory  and  time.  An  alternative  procedure 
would  be  to  call  NPCHEMIN  (and  NPREPEQ  if  necessary)  once  for  each  reaction  set,  and 
to  write  the  variables  stored  in  the  common  blocks  to  separate  Hies  or  logical  records.  Then 
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on  crossing  the  boundary  between  flow  regions,  a  single  read  statement  would  provide  ail 
the  reaction  rates  and  properly  indexed  material  properties  relevant  to  the  domain.  The  user 
would  still  need  to  develop  his  own  method  to  properly  pass  information  from  domain  to 
domain  as  NEQPAK  does  not  provide  such  a  capability. 

Table  6.  Calling  Sequence  for  Subroutine  NPREPEQ 
SUBROUTZMB  llI>ItSPEq<FSCOMP  .RHDFS . T1  ,T2 , PI , P2. HPRS  .ITHP .  lEq) 

Outputs  ths  squilibriuB  conposition  at  equally  spaced  temperatures  and  pressures. 

Sutooutine  VPCHEKIH  must  be  called  first  to  set  up  species  property  arrays. 

IIPUT 

FSCOMP(IIXBS)  —  Free-streau  composition  (kg-moles/n^) 

REOFS  —  Free-stream  density  (kg/m^) 

T1,T2  —  Minimum  and  nazumim  (K) 

temperature 

Pt,P2  —  Minimum  and  mazinum  pressure  (Pa) 

WPRS  —  lumber  of  pressures  — 

HTNF  —  lumber  of  temperatures  — 

lEQ  —  Fortran  unit  number  for  — 

output 

KQ7E 

T1  may  be  greater  or  less  than  T2.  PI  may  be  greater  or  less  than  P2. 

3.3  CHEMICAL  SOURCE  TERMS 

The  species  source  terms  vv^  in  Eq.  (3)  and  their  derivatives  with  respect  to  the  independent 
variables  7^  and  are  evaluated  in  Subroutine  NPSOURCE.  The  gain  and  loss  terms, 
andC.r,f,  defined  in  Sec.  2.4  are  returned  in  addition  to  iV^,  since  some  specialized  integration 
algorithms  treat  the  two  terms  differently.  The  calling  sequence  is  shown  in  Table  7. 
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Table  7.  Calling  Sequence  for  Subroutine  NPSOURCE 
SUBROUTIVE  IPSOURCEtGAM, TEMP, DERI V , LOW, IP , IlC, TMII , GAZI . LOSS . VOGT , DHDTM ,DUDG) 


IIPOT 


GAN(MXIS,IP) 

7. 

Concentration  of  species  s 
at  point  n 

(kg-noles/n^) 

TEMP(rT*HP) 

T, 

Tenperatnre  of  type  q  at 
point  n 

(K) 

DERZV 

Flag  to  control  calculation 
of  derivatives 

— 

LQV.IP.IIC 

— 

Do-loop  controls 

— 

TMII 

Tmin 

Niniann  tes^eratnxe  for 
calculating  reactions 

OUTPUT 

(K) 

6AII(KXlR,HXaS.IP) 

Gr,t 

Incrsaso  of  species  s  at  n 
by  reaction  r 

(kg-noles/s/n^) 

L0SS(NX1R.MXIS,IP} 

Decrease  of  species  •  at  n 
by  reaction  r 

(kg-Boles/ b/b^ ) 

tfDOT(MXlS,IP) 

let  change  of  species  s  at 
point  n 

(kg-aoles /s/b^ ) 

DVDTN(NXIT,HXI5,IP) 

dWtfSTM 

Derivative  of  UDOT  srt. 
tenperature  m  at  n 

(kg-noles/s/K/B^ ) 

DUDG(MXIS,MXIS,IP) 

dwjdtj 

Derivative  of  VDQT  srt. 
concentration  of  species  j 

at  n 

(■-^) 

3.3.1  The  TEMP  Array 

Both  jC,.  j  and  ,  involve  the  forward  rate  coefficients  defined  by  Eq.  (28)  in  terms 
of  a  “generic”  temperature  which  may  differ  from  reaction  to  reaction  and  may  well 
be  some  user-defined  function  as  discussed  in  Sec.  2.4  above.  These  chemical  temperatures 
are  supplied  to  NPSOURCE  through  the  TEMP  array. 

TEMP  is  a  singly-dimensioned  array  of  length  NT  times  NP  where  NT  is  the  number 
of  different  temperatures  used  in  a  particular  problem  and  NP  is  the  number  of  grid  points. 
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The  variable  NT  is  determined  by  NPCHEMIN  from  the  input  reaction  data  and  stored 
in  common  block  NEQCOM.  The  value  of  NT  is  returned  to  the  flow  program  via  the 
argument  NTMP. 

NPSOURCE  evaluates  the  forward  rate  for  reaction  number  IR  at  grid  point  IP  using 
the  temperature  TEMP(IT),  where 

IT  ■  (TTYPEF(IR)-1)+NP  +  IP 

and  TTYPEF(IR)  and  TTYPER(IR)  are  arrays  used  to  store  the  information  on  which  generic 
temperature  is  used  to  evaluate  the  forward  and  backward  rates.  TTYPEF  and  TTYPER 
are  stored  in  the  common  block  /NEQCOM/  (Appendix  C)  and  correspond  to  the  variables 
TYPEF  and  TYPER  read  by  NPCHEMIN  from  the  reaction  data  file  (Appendix  F).  By 
convention,  TYPE  =  1  is  assigned  to  the  translational  temperature  of  the  heavy  particles 
and  TYPE  =  2  is  assigned  to  the  translational  temperature  of  the  free  electron.  Thus,  NT 
is  greater  than  or  equal  to  2  except  for  a  nonionized  gas  in  thermal  equilibrium.  For  the 
two-temperature  model,  the  vibrational  temperature  and  the  electron  temperature  are  equal 
and  are  both  assigned  TYPE  =  2. 

The  TEMP  array  is  also  used  by  Subroutine  NPSVC  (Sec.  3.4.4)  when  the  user  models 
vibration-chemistry-vibration  coupling  with  one  of  the  models  described  in  Secs.  2.4.2  and 
2.5.4.  When  this  is  the  case,  if  Model  3  were  chosen,  then  TYPE  =  3  would  correspond 
to  r„.  If  Model  4  were  chosen,  then  TYPE  =  3  through  TYPE  =  1  +  would  correspond 
to  the  vibrational  temperatures  of  the  species  which  the  user  specified  to  be  in  thermal 
nonequilibrium.  The  temperatures  must  be  entered  into  the  TEMP  array  in  the  order  in  which 
the  vibrators  appear  in  the  FSCOMP  array  on  return  from  NPCHEMIN  (Sec.  3.2.1  ).  The 
user  may  wish  to  write  a  preprocessor  code  which  calls  NPCHEMIN  to  determine  this  order. 

Any  user-defined  temperatures,  such  as  that  suggested  by  Park  (Ref.  13)  (Sec.  2.4),  are 
the  final  entries  into  the  TEMP  array  and  would  be  assigned  TYPE  >  3. 

3.3.2  Rate  Coefficients 

The  forward  reaction  rate  coefficient  is  modeled  with  the  modified  Arrhenius  form,  Eq. 
(28),  namely 

kf  =  ArTf^exv{-Cr(T,) 


and  may  be  any  one  of  the  temperatures  defined  by  the  various  kinetic  models. 
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3.3.3  Calcnlation  of  Reverse  Rates 

The  method  of  calculating  the  reverse  rate  for  the  rth  reaction  is  communicated  to 
NPSOURCE  through  the  flag  REV(r)  which  has  a  value  of  0,  1,  or  2  as  determined  by 
Subroutine  NPCHEMIN  from  the  reaction  data  file  (see  Appendix  E).  Unless  otherwise 
specified,  REV(r)  =  1  and  Eqs.  (29)  and  (30)  are  applied  using  a  calculated  in  Subroutine 
NPGIBBS  from  the  polynomial  fits  of  Eq.  (30). 

If  REV(r)  =  2  ,  the  reverse  rate  coefficient  is  calculated  directly  using  the  curve  fits 

DrTf’cxpi-FrIT,)  (67) 

which  are  entered  for  each  reaction.  The  reverse  rate  is  evaluated  using  the  temperature  type 
TTYPER(r). 

If  REV(r)  =  0 ,  the  reverse  rate  is  not  evaluated  at  all.  This  is  useful  in  modeling  irreversible 
reactions  such  as  those  found  in  certain  global  combustion  models. 

3.3.4  Source  Term  Derivatives 

Since  the  source  terms  in  Eqs.  (3)  and  (6)  tend  to  be  “stiff,”  the  use  of  implicit  methods 
for  solving  the  conservation  equations  has  become  common  practice.  However,  many 
algorithms  for  integrating  “stiff”  systems  of  differential  equations  require  the  evaluation 
of  Jacobian  derivatives.  Since  the  evaluation  of  these  elements  is  time  consuming,  it  is 
sometimes  advantageous  to  omit  the  re-evaluation  of  a  Jacobian  for  several  iteration  cycles 
or  marching  steps.  A  common  approximation  in  many  algorithms  is  to  use  only  the  diagonal 
terms  of  the  Jacobian,  i.e., 

The  variety  of  possible  thermodynamic  models  and  conservation  variables  makes  it 
impractical  to  furnish  the  source  term  derivative  values  in  exactly  the  form  required  by  each 
solution  algorithm.  Instead,  NPSOURCE  and  all  other  NEQPAK  source  modules  calculate 
the  source  term  derivatives  with  respect  to  the  independent  variables  supplied  to  the  module. 
If  the  logical  variable  DERIV  =  .TRUE. ,  NPSOURCE  returns  the  derivatives  dw/dyj  and 
div/dTjj,  where may  be  any  one  of  the  temperatures  in  the  TEMP  array. 

Users  must  then  apply  the  chain  rule  to  form  the  Jacobian  elements  to  suit  their  own  needs. 
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The  derivatives  which  enter  the  Jacobian  are  not  entirely  independent.  In  the  first  place, 
conservation  of  atoms  imposes  linear  constraints 

JV, 

^  ^  ®t,m Wtii  —  0.0,  I  —  1,.  • . ,  NfUg 

m=l 

where  is  the  number  of  atoms  of  element  /  in  species  m,  and  is  the  number  of 
different  atom  types.  Second,  conservation  of  mass  requires  that 


(68) 

where  q  is  the  ordinary  global  density  and  is  the  molecular  weight  of  species  s.  Finally, 
if  mass  fractions  or  mole  fractions  X,  =  7j/Ej  are  used  in  formulating 

the  conservation  equations,  then  these  fractions  must  sum  to  1 .0.  References  14  and  47  discuss 
the  wisdom  of  using  the  foregoing  constraints  to  replace  one  or  more  of  the  species  conservation 
equations.  Since  there  are  so  many  choices  available,  NPSOURCE  evaluates  the  derivatives 
as  if  these  constraints  did  not  exist.  However,  if  the  reactions  are  input  correctly,  both 
and  its  derivatives  will  automatically  satisfy  the  constraints.  Thus,  users  may  safely  manipulate 
the  values  returned  by  NPSOURCE  to  suit  their  own  integration  strategy. 

3.3.5  Vibrational  Coupling 

If  any  reaction  is  input  with  a  non-zero  value  for  the  flag  CTYPE(r) ,  NPCHEMIN  sets 
the  logical  variable  COUPLE  =  .TRUE,  and  NPSOURCE  calls  Subroutine  NPCOUPLE 
to  calculate  a  vibrational  coupling  factor  for  the  reaction.  As  discussed  in  Sec.  2.4.2,  this 
method  of  modeling  the  effect  of  vibrational  excitation  on  chemical  reaction  rates  assumes 
that  the  quoted  rate  coefficients  were  measured  for  a  thermal  equilibrium  situation  and 
multiplies  Eq.  (28)  by  a  function  of  translational  and  vibrational  temperatures.  NEQPAK 
extends  this  model  slightly  to  allow  separate  coupling  factors  iy  and  for  the  forward  and 
reverse  directions.  These  coupling  factors  are  calculated  in  Subroutine  NPCOUPLE. 

Coupling  models  proposed  by  Olynick  and  Hassan  (Ref.  11),  by  Landrum  and  Candler 
(Ref.  12),  and  by  Knab,  Friihauf,  and  Jonas  (Ref.  10)  are  implemented  in  NPCOUPLE. 
NPCOUPLE  uses  the  array  CTYPE(r)  to  determine  which  coupling  model  to  apply  to  reaction 
number  "r."  CTYPE(r)  corresponds  to  the  variable  KPL  read  by  NPCHEMIN  from  the 
reaction  data  file  (assigned  to  unit  number  INPR).  In  the  reaction  data  Hie,  the  user  should 
set  KPL  =  1  for  Olynick-Hassan  coupling,  KPL  =  2  to  select  the  Landrum-Candler  model, 
or  KPL  =  3,  4,  or  S  to  choose  one  of  the  coupling  models  developed  by  Knab,  Friihauf, 
and  Jonas.  The  user  should  set  KPL  =  0  when  using  the  one-temperature  model  or  when 
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modeling  the  effects  of  thermal  nonequilibrium  on  reaction  rates  with  Park’s  model.  The 
other  parameters  needed  by  a  particular  coupling  model  are  entered  via  the  G,  H,  and  U 
input  variables.  During  the  reaction  processing  performed  by  NPCHEMIN,  the  values  of 
G  »  H,  and  U  are  loaded  into  the  common  arrays  GK(t),  HK(r),  and  UK(r). 

CTYPE  =  1:  Olynick — Hassan  Coupling  The  Olynick-Hassan  coupling  model  was  developed 
for  dissociation  reactions  of  the  form 

where  molecule  AB  is  a  rotationless  harmonic  oscillator  and  Af  is  an  arbitrary  third  body. 
In  their  analysis,  Olynick  and  Hassan  assumed  that  the  vibrational  energy  of  AB  is  given 
by  a  Hinshelwood  distribution.  Since  the  Hinshelwood  distribution  does  not  allow  for 
population  inversions,  the  Olynick-Hassan  model  is  not  valid  for  expanding  flows. 

Let  z  =  Cv/fTZT’y)  and  define 

x  =  {l~TJT)Ti,JT^ 

where  7^^  is  the  dissociation  temperature  of  AB.  Then,  the  forward  coupling  factor  is 

_  »T(z,z) 

where  r(2,x)  is  the  incomplete  gamma  function.  The  reverse  coupling  factor  is 

9r  =  1.0 


The  incomplete  gamma  function  is  evaluated  by  the  function  NPSGAM  using  an  algorithm 
from  Press,  et  al.  (Ref.  27).  Input  of  the  Olynick-Hassan  model  into  the  reaction  mechanism 
file  is  illustrated  by  the  reaction  set  “hassan”  in  Appendix  F-1. 

CTYPE  =  2:  Landrurn'Candler  The  Landrum-Candler  model  was  developed  from  an  intensive 
numerical  investigation  of  the  dissociation  of  N2.  Its  applicability  to  other  diatomic  species 
has  not  been  investigated,  so  the  generalization  coded  here  may  not  be  warranted.  The  essence 
of  the  model  is  to  replace  the  dissociation  of  the  diatomic  molecule  AB  by  the  two-step  reaction 

AB  +  M  ^  AB*  +  M 
AB*  +  M  ^  A  +  B+M. 

In  the  above  reactions,  AB*  represents  ihtAB  molecules  with  vibrational  quantum  numbers 
V  where  '^d.AB  dissociation  temperature  of  AB,  '^d.AB 

=  Based  on  the  results  of  their  numerical  investigation,  Landrum  and  Candler 
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assumed  that  the  virtual  species  AB*  is  in  thermal  equilibrium,  but  AB  has  a  vibrational 
temperature  Ty  ^  T.  Landrum  and  Candler  argued  that 


for  the  reaction 


=  exp(/3dxfl 


T  -  rt,.exp[-(i>e  -  \)9ABfTv\ 

T„r  Z{vc6AB,Ty) 


>1J3  +  M  ^  A5*  +  M 


In  the  definition  of  ♦y  above,  is  the  characteristic  vibrational  temperature  of  AB,  and 
Z  is  the  vibrational  partition  function  for  >45. 2is  computed  under  the  assumption  that  the 
vibrational  energy  mode  of  >45  is  described  by  a  harmonic  oscillator  truncated  at  quantum 
number  -  1.  Zis  calculated  by  Subroutine  NPART.  The  parameter  5  =  GK(IR)  is 
supposed  to  correct  for  the  anharmonicity  of  >45.  Landrum  and  Candler  were  able  to  obtain 
impressive  agreement  with  their  full  simulation  results  for  N2  by  taking  iS  =  S.O  and  increasing 
the  Millikan-White  (Ref.  28)  vibrational  relaxation  time  t  [Eq.  (46)]  by  a  factor  of  25. 


The  reactions  in  Eq.  (33)  are  entered  into  the  reaction  file  in  the  usual  way  as  illustrated 
by  the  reaction  set  “candler”  shown  in  Appendix  F-1.  In  the  last  reaction  of  this  set,  the 
electron  temperature  (TYPER  =  2)  is  assumed  to  control  the  recombination  of  NO  +  and  E  - , 


NPCHEMIN  assigns  the  same  thermodynamic  and  transport  properties  to  the  virtual 
molecule  AB*  as  to  AB  so  it  is  not  necessary  for  the  user  to  make  any  special  entries  in 
the  database.* 


CTYPE  3,4,5:  KNAB,  FRUHAUF,  and  JONAS  Knab,  Friihauf,  and  Jonas  (Ref.  10) 
have  extended  the  concept  of  vibrational  coupling  to  reaction  types  other  than  dissociation. 
Table  8  shows  for  each  type  of  reaction  the  coupling  factors  ’F,(r,7\,)  recommended  by  Knab, 
et  al.  and  the  CTYPE  values  assigned  to  the  reactions  by  NEQPAK. 

The  input  parameter  =  HK(r)  +  UK(r)  •  T  measures  the  probability  that  a  molecule 
in  vibrational  state  /  will  dissociate.  If  0  =  oo,  dissociation  is  equally  probable  from  all 
vibrational  levels.  The  values  for  HK  and  UK  are  read  from  the  reaction  mechanism  file. 
Knab,  et  al.  also  introduce  an  effective  activation  temperature  =  aC,.  for  reaction  r.  An 
appropriate  value  for  a  is  supplied  via  the  input  variable  G.  If  one  sets  a  =  1 ,  the  formulation 
of  Eq.  (37)  reduces  to  that  of  Treanor  and  Marrone  (Ref.  25). 

Reaction  set  “fruhauf”  in  Appendix  F-1  illustrates  the  input  for  all  three  reaction  types 
listed  in  Table  8. 


*  In  order  to  completely  implement  Landrum  and  Candler’s  method,  the  user  would  have  to  redo 
the  curve  fits  for  both  AB  and  AB*  using  the  thermal  models  of  the  molecules  defined  in  Ref.  12. 
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Table  8.  Coupling  Factors  for  Various  Reaction  Types 


CTYPE 

Reaction 

Forward  Rate 

Reverse  Rate 

3 

AB  +  M^A  +  B  +  M 

K 

4 

AB  +  C  ^  AC  +  B 

5 

A  +  B  ’^AB-  + 

kf(T) 

3.4  VIBRATIONAL  ENERGY  SOURCE  TERMS 

The  species  vibrational  energy  is  affected  by  inelastic  collisions  with  heavy  molecules  and 
with  electrons.  During  a  heavy  particle  collision  there  may  be  an  exchange  of  vibrational 
energy.  The  species  vibrational  energy  is  also  affected  by  an  increase  or  decrease  in  the  number 
of  vibrationally  excited  molecules  through  chemical  reaction.  These  interactions  are  modeled 
by  the  Subroutines  NPEVDOT,  NPSVE,  NPSVV,  and  NPSVC,  respectively,  as  discussed 
below.  Interaction  with  the  electromagnetic  field  is  not  modeled  at  the  present  time.  All  inputs 
into  the  vibrational  energy  source  term  subroutines  dealing  with  species-dependent  Quantities 
are  dimensioned  with  the  parameter  MXNS  and  are  indexed  according  to  species  number 
in  the  order  returned  by  NPCHEMIN.  On  output,  all  arrays  dealing  with  the  vibrational 
energy  source  term  and  its  derivatives  are  dimensioned  with  the  parameter  MXNV  and  are 
indexed  according  to  vibrator  number  in  the  order  in  which  the  vibrators  appear  in  the 
FSCOMP  array  on  return  from  NPCHEMIN.  If  the  logical  variable  DERIV  =  .TRUE., 
Subroutine  NPSVC  will  return  derivatives  of  w,  with  respect  to  7^,  T,  and  The  other 
source  modules  will  return  appropriate  derivatives  regardless  of  DERIV.  If  Model  2  or  Model 
3  is  chosen,  Ejl'',  7,+|^fgtI>,  and  Eflvi  7;^.^g3u/9r„  are  returned  in  the  positions  for  vibrator 
number  1  where  is  the  number  of  elements  (including  free  electrons)  present  in  the 
mixture.  All  the  vibrational  energy  source  term  routines  return  the  time  rates  of  change  of 
energy  densities  in  terms  of  IV/ m^. 

3.4.1  V-T  Interaction 

Landau  and  Teller  (Ref.  29)  showed  that  if  the  vibrational  energy  levels  of  a  diatomic 
atom  are  equally  spaced  (harmonic  oscillator  assumption)  and  that  if  an  inelastic  collision 
with  another  atom  induces  only  single-level  transitions  among  the  vibrational  energy  levels 
then  the  time  rate  of  change  of  the  vibrational  energy  as  it  relaxes  towards  equilibrium  is 
given  by 

51,-r  ^ 
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where  the  ef  is  the  vibrational  energy  computed  with  the  translational  temperature  and  t, 
is  the  so-called  relaxation  time.  SJ'^is  computed  in  NEQPAK  in  Subroutine  NPEVDOT. 
The  calling  sequence  for  NPEVDOT  is  shown  in  Table  9.  On  input,  NPEVDOT  requires 
the  species  vibrational  energies  and  specific  heats  at  constant  volume  When  calling 
NPEVDOT  with  Model  2,  the  user  must  first  call  Subroutine  NPTHRM  (Sec.  ^6)  with  MODE 
=  3  and  Tg  =  to  separate  the  contributions  of  the  electron-electronic  excitation  energy 
modes  from  the  vibrational  energy  modes.  For  all  species,  ej  =  AJ  and  CJ Strictly  speaking, 
Eq.  (42)  is  valid  only  for  a  harmonic  oscillator,  but  it  is  applied  to  anharmonic  oscillators 
in  Subroutine  NPEVDOT  by  calling  NPTHRM  to  calculate  the  equilibrium  vibrational 
energy  e^*. 

The  species  relaxation  time  Tj  is  evaluated  in  Subroutine  NPRELAX  which  is  called 
from  NPEVDOT. 

Subroutine  NPLAXIN,  which  is  called  by  NPCHEMIN,  reads  the  necessary  relaxation 
data  from  unit  INREL.*  If  NPLAXIN  does  not  find  data  for  a  particular  interaction  pair, 
the  subroutine  generates  the  necessary  coefficients  by  making  use  of  the  correlation  of  Milliken 
and  White,  Eq.  (46). 

3.4.2  V-E  Interaction 

An  encounter  between  a  free  electron  and  a  molecule  may  result  in  the  transfer  of  kinetic 
energy  from  the  electron  to  a  vibrational  mode  of  the  heavy  particle.  The  probability  for 
such  an  interaction  is  an  extremely  irregular  function  of  the  electron  kinetic  energy. 
Nevertheless,  Lee  (Refs.  8  and  30)  has  been  able  to  model  the  process  with  an  equation  of 
the  Landau-Teller  type  [Eq.  (42)] 


where  is  the  vibrational  energy  evaluated  at  the  electron  temperature  and  the  relaxation 
time,  Tgg,  is  given  by  ^ 

"  n^il-e-o^f^-yfko,vv^dv 


*  Provision  has  been  made  to  allow  NPLAXIN  to  read  this  block  from  a  different  unit  from  that 
used  for  the  rest  of  the  thermal  properties  database  to  allow  some  experimentation  with  relaxation 
rates  without  the  risk  of  corrupting  these  data. 
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Table  9.  Calling  Sequence  for  Subroutine  NPEVDOT 
SUBROUTINE  ipevdotCtt.te.esv.cvsv.gan.p.loh.ip.ivc.devdt.evdotdt.evdotdg.evdqttv) 


INPUT 

TT(HP) 

T 

Transla.tloiial  taup«ratiir«  at 
point  n 

CK) 

tecbp) 

Blactronic  taDpsratnra  at  u 

(K) 

ESV(NXBS,1P) 

< 

Vibrational  energy  oT 
species  s  at  n 

(j/kg-BOlo) 

CV&V(KXIS,NP) 

Vibrational  specific  heat  of 
apacles  s  at  n 

(j/kg-nole/K) 

GAKCMINS.IP) 

7» 

Concentration  of  species  s 
at  point  n 

(kg-noles/n^) 

PfSP) 

P 

Static  pressure  of  the 
nixtnre  at  n 

(Pa) 

LOW, NP, INC 

Do- loop  controls 

FRQX  COHHON 

MODEL 

2,3,4 

Themal  nodel  definition 

OUTPUT 

■  " 

DEVDT(HXNV,NP) 

Tine  rate  of  change  of 
vibrational  energy  of 
nixtnre  or  Individual 
species  at  n 

(W/b3) 

EVDaTOT(NXHV,NP) 

dS“-'^/dT 

Derivative  of  DEVDT  vrt. 
translational  tenperature  at 

(¥/nVK) 

EVDOTDG(MXirV,NP) 

dsr^ldyi 

JR 

Derivative  of  DEVDT  vrt. 
the  concentration  of  species 
j  at  n 

(V/kg-Bole) 

EVDQTTV(HXNV,VP) 

dsi-'^ldT,,. 

Derivative  of  DEVDT  art. 
vibrational  tenperature  at  n 

(W/n^/K) 

where  is  the  electron  number  density,  is  the  characteristic  vibrational  temperature  of 
species  s,  and  ^  is  the  rate  constant  for  a  transition  from  the  vibrational  ground  state  to 
vibrational  level  v. 

The  terms  are  calculated  in  Subroutine  NPSVE.  The  calling  sequence  for  NPSVE 
is  shown  in  Table  10.  Subroutine  NPSVE  uses  polynomial  fits  to  the  natural  logarithm  of 
the  integral  term  of  Eq.  (49)  to  calculate  the  relaxation  rate  for  electron  impact.  The  rate 
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data  from  Lee  (Ref.  30),  Huo,  et  al.  (Ref.  48),  Kieffer  (Ref.  49),  and  Stinker  and  Ali 
(Ref.  50)  were  integrated  numerically  and  fitted  at  AEDC. 

Table  10.  Calling  Sequence  for  Subroutine  NPSVE 
SUBROUTIVB  IPSVE (TE . G&M , CVSV , ESV ,LDU , IP , I¥C . SVE , SVETE , SVEG . SVFTV ) 

IIPOT 


TE(IP) 

Electronic  tanpsratnrc  at  n 

(K> 

GAH(NUrS.BP) 

7« 

Concantration  ol  apaciaa  s 
at  point  n 

(kg-molea/m^) 

ESV(HXIS,BP) 

Vibrational  anargy  ol 
apaciaa  a  at  n 

(j/kg-nola) 

CVSV(HXIS,IP) 

Vibrational  spacific  heat  of 
apaciaa  a  at  n 

(j/kg-nole/K) 

LOU. BP, INC 

Do-loop  controla 

FROM  COMMON 

EVCDEF(£.HXBS) 

— 

Coafficianta  In 
croBB-aaction  polynomial 

— 

MODEL 

2,3.4 

Thamal  nodal  daf  inition 

OUTPUT 

SVB(MXBV,IP) 

Excitation  rata  by  electron 
impact  at  n 

(u/m^) 

5VETE(MXIV,BP} 

dS)-^ldT^ 

Darivativa  of  SVE  art. 
elactronlc  taaparature  at  n 

(W/m^/K) 

SVEGCMXNV.NP) 

Darivativa  of  SVE  art,  tba 
concantration  of  apaciaa  j 
at  n 

(V/kg-mola) 

SVETV(MXBV,IP) 

Darivativa  of  SVE  art. 
vilnratlonal  tamparatnra  at  n 

(W/m®/K) 

Since  the  vibrational  energy  gained  by  the  molecule  is  lost  by  the  electron,  NPSVE  should 
not  be  called  for  Model  2  in  which  vibrational  and  electron  energy  are  pooled. 


3.4.3  V-V  Interaction 

It  is  possible  for  two  molecules  to  exchange  vibrational  energy  during  a  collision,  but 
this  process  is  unlikely  unless  the  molecules  happen  to  have  vibrational  levels  Vj  and  Vj  with 
nearly  the  same  energy.  This  vibrational  exchange  may  be  considered  to  be  a  reaction  such  as 
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E2(l)  +  ^20(000)  =  ff2(0)  +  ^20(001) 


which  describes  the  InteTchange  of  vibrational  energy  between  H2  and  one  of  the  modes  of 
H20.  Rate  data  for  such  V-V  reactions  are  available  (Ref.  31),  but  considering  level-to-level 
transitions  of  the  kind  indicated  is  computationally  expensive.  Therefore,  NEQPAK  has 
adopted  the  model  of  Candler  and  MacCormack  (Ref.  15)  in  the  form  proposed  by  Knab, 
et  al.  (Ref.  10).  The  crude  model  described  by  Eq.  (50)  is  evaluated  in  Subroutine  NPSVV 
with  the  calling  sequence  of  Table  11. 

3.4.4  V-C  Interactions 


It  is  clear  that  the  dissociation  of  a  vibrating  molecule  removes  its  vibrational  energy  from 
the  species  pool  and  that  the  formation  of  a  new  molecule  must  add  some  vibrational  energy 
to  the  pool.  NEQPAK  offers  the  user  several  models  of  vibration-chemistry  coupling  to  choose 
from.  The  models  are  coded  in  Subroutine  NPSVC  whose  calling  sequence  is  shown  in  Table 
12.  Calls  to  NPSVC  must  be  preceded  by  calls  to  NPRELAX  (Table  13)  and  NPSOURCE. 
NPSOURCE  must  be  called  to  obtain  values  for  Wj,  £,  j,,  and  dw^/dlj,  and  dw^/dT^. 
The  call  to  NPRELAX  is  necessary  to  get  values  for  t,,  dr/dT,  and  dr/dyj. 


Non-preferential  Dissociation  The  simplest  assumption  is  that  the  energy  gained  or  lost  during 
the  reaction  is  just  the  average  vibrational  energy  of  the  species,  i.e., 

5,”-'  = 

This  non-preferential  dissociation  assumption  is  appropriate  for  Model  2  and  also  for  Olynick- 
Hassan  coupling  (Ref.  1 1)(CTYPE  =  1).* 


Landrum-Candier  Coupling  The  Landrum-Candier  theory  (Sec.  2.4.2)  assumes  that  vibrational 
levels  V  <  Vp  relax  according  to  the  Landau-Teller  theory  but  the  upper  vibrational  levels 
(the  virtual  species  AB*)  are  in  thermal  equilibrium-  The  critical  level  is 


Vc 


Mg  -  T 
Oab 


+  1 


•  Since  all  the  terms  needed  to  compute  the  term  defined  in  Eq.  (51)  are  available  by  calling 
NPSOURCE  (Sec.  3.3)  and  NPTHERM  (Sec.  3.6),  users  might  choose  to  construct  5^“*’  and  the 
appropriate  derivatives  within  the  flow  solver  as  an  alternative  to  calling  NPSVC. 
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where  is  the  dissociation  temperature  of  AB  and  is  the  characteristic  vibrational 
temperature  for  species  AB. 

Only  defined  by  Eq.  (53)  and  its  derivatives  with  respect  to  temperature, 
concentration,  and  vibrational  energy  are  returned  by  Subroutine  NPSVC.  The  Landau-Teller 
term  at  the  beginning  of  Eq.  (52)  is  calculated  separately  by  a  call  to  NPEVDOT  (Table  9). 

Table  11.  Calling  Sequence  for  Subroutine  NPSVV 
SUBROUTIIE  HPSVVCTT ,  GIK , EVS  ,CVSV  ,L0W , HP , IIC .SVV.SVVT,  SWG , SWTV) 

IIPUT 


TT<KP) 

T 

Tjranslatlonal  tanparature  at 
point  n 

(K) 

GAMCRZIS.HP) 

7i 

Concentration  ot  epeciaa  a 

at  n 

(kg-noles/m^) 

EVS(iaiS,VP) 

< 

Vibrational  anargy  ot 
spacias  a  at  n 

(j/kg-nole) 

CVSV(MXHS,HP) 

Vibrational  spaeifie  heat  of 
spacias  s  at  a 

(J/hg-mola) 

L0V,HP,I1C 

Do-loop  controls 

FROM  CDHItOI 

THET(HXHS) 

Characteristic  vibrational 
tenperatnre  of  species  a 

(K) 

SLJ(1[XIS,HXHS) 

^i,j 

Lanard- Jones  collision 

diaiseter 

OUTPUT 

(m) 

SVV(MXHV,HP) 

gv-v 

Vibrational  energy  source 
tana  at  n 

(W/n®) 

SVVTCHMV.HP) 

Derivative  of  SW  art. 
translational  tenparature  at 

n 

<W/ii®/K) 

SVVG<HXHV,HX¥S,HP) 

dsrvdii 

Derivative  of  SW  srt. 
concentration  of  species  j 

(W/kg-»ole) 

SWTVdfXIV.HXIV  ,XP) 

Derivative  of  SW  nt. 
vibrational  tenparatnre  of 
species  j 

(V/ia®/K) 

KFJ  Coupling  The  coupling  theory  of  Knab,  et  al.  (Ref.  10)  is  invoked  for  reactions  with 
CTYPE  =3,4,5.  Unlike  the  other  coupling  theories  described  above,  in  the  KFJ  theory  (Ref. 
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10)  the  vibrational  coupling  factor  depends  on  both  the  translational  and  the  species  vibrational 
temperatures.  Thus,  NPSVC  returns  5]'“'^  along  with  its  derivatives  with  respect  to 
translational  temperature,  species  concentrations,  and  species  vibrational  temperatures  when 
KFJ  coupling  is  invoked. 

Table  12.  Calling  Sequence  for  Subroullne  NPSVC 

IPSVC(TBKP.GAK,G&ZI,LOSS,UDOT,EVS,CTSV,TAU,DVDG,DWl}TN,DTAUT.I>TiU6,DEaZT,LOU,IP.IIC,SVC,SV(rr,SVCQ) 

ZIPtJT 


TBIIP(IT*BP) 

Taiperature  of  tjpo  q  at 
point  a 

(K) 

GAMCntlS.lP) 

7* 

Cone oatrat ion  of  cpociaa  ■ 
at  point  a 

(kg-Bolas/n^) 

GAlKMXIR.HZIS.IP) 

Incraaao  of  apoeios  a  at  a 
bp  raactlon  r 

(kg-Bolaa/a/n^) 

LOSSCBXni.lIltBS.IP) 

Oacraaae  of  apacias  a  at  a 
bj  raactlon  r 

(kj^^las/a/n^) 

VD0TCN11S,IP) 

«?* 

lat  changa  of  apaciaa  a  at 
point  n 

(kg-«olaa/a/«®) 

BVS(MUS,I7) 

Vitaational  anargp  of 

apaciaa  a 

CVSV(HXIS.IP) 

Vibrational  boat  capacity  of 

apacias  a 

(j/Cn^E)) 

L0V,iP,irc 

— 

Do-loop  controls 

— 

TAOCMXIS.HP) 

r. 

Raleucatlon  time  of  apaciaa  a 

at  n 

(a) 

DTAUTfXXVS.IP) 

drJdT 

Darivativa  of  TAU  aith 
raspect  to  tranalatioBal 
tai^aratiira  at  a 

(a/I) 

OTAUG(MnS.I(XIS,IP) 

dr^fd’^j 

DariTatira  of  TAU  with 
raspect  to  coneantration  of 
apaciaa  j  at  n 

(an^/kg-nola) 

DWDG(NZlS,iaiS,IP) 

Dari vat iva  of  HDOT  vlth 
respect  to  coneantration  of 
apacias  j  at  n 

(a-1) 

D«mi(Hxrr,Nzis,ip) 

dwgf&Tm 

Derivative  of  HDOT  vitb 
raapact  to  taaparatura  of 
type  ■ 

(kg-nolaa/ s/ib^/K  ) 
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Table  12,  Calling  Sequence  for  Subroutine  NPSVC  (Concluded) 

FRCm  COMHOI 


TEET(MXI5) 

Characteristic  vibrational 
teoperatore 

(K) 

DISOC(HXItS} 

A. 

Dissociation  energy  of 
species  s 

(j/kg-mole) 

CK(KXXR) 

Cr 

Activation  energy  of 
reaction  r 

(K) 

GK.HK,IIK(NXXK) 

— 

Paraneters  for  reaction  r 

— 

CTYPECHXHR) 

— 

Coupling  type  for  reaction  r 

— 

NC 

— 

NuIbud  value  of  ctype 

— 

NT 

— 

Hnnber  of  temperature  types 

OITTPDT 

— 

DERIV 

— 

Flag  to  control  calculation 
of  derivatives 

— 

SVC(MXNV,NP) 

Rate  of  increase  of 
vibrational  energy 

svcKmxnv.wp) 

dSr^fdT^ 

Derivative  of  SVC  vrt 
teBperature  m 

(W/n®/K) 

SVCGCHXIV.MXIS.NP) 

ds--/djj 

Derivative  of  SVC  vrt 
concentration  of  species  j 

CV/kg-Dole) 

46 


AEDC-TR-93-20 


Table  13.  Calling  Sequence  for  Subroutine  NPRELAX 

SUBHOITTIIE  NPRELAKTT , GAM, P, LOW. IP , IIC ,TAU,DTAUT .DTiUG) 


lIPBT 

TT(SP) 

T 

Translational  temperature  at 
point  n 

(K) 

GAM(HXIS.IP) 

7* 

Concentration  ol  species  s 
at  point  n 

(kg-noles/m^) 

P(»P) 

P 

Pressure  at  n 

(Pa) 

LOW.IP,IIC 

— 

Do -loop  controls 

OUTPUT 

— 

TAU(MXBS,1P) 

Relaxation  tine  of  species  a 
at  n 

(s) 

DTAOTCHXHS.HP) 

drJdT 

Derivative  of  TAU  vrt  to 
translational  temperature  at 

et 

(s/K) 

DTAUG(NnrS,NXIS,IP) 

4A 

Derivative  of  TAU  vith 
respect  to  concentration  of 
species  j  at  n 

(a-n®/kg-mol* 
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3.5  ELECTRON-ELECTRONIC  ENERGY  SOURCE  TERMS 

The  energy  of  the  free  electron  is  affected  by  elastic  and  inelastic  collisions  with  heavy 
particles,  by  electron  production,  by  interaction  with  applied  and  induced  electric  fields,  and 
by  radiation  processes.  The  energy  of  electronic  excitation  may  be  lost  to  radiation.  Of  these 
effects,  only  those  dealing  with  collisions  and  electron  production  through  chemical  reaction 
are  included  with  the  NEQPAK  source  terms  [ij'  of  Eq.  (6)].  Each  of  the  source  terms 
discussed  in  Sec.  2.6  is  calculated  in  a  separate  module,  but  a  single  call  to  Subroutine 
NPEEDOT  (Table  14)  will  return  w*  defined  in  Eq.  (9).  NPEEDOT  also  returns  the 
derivatives  of  with  respect  to  translational  and  electron-electronic  temperatures  in  the  array 
EEDOTTM.  The  derivatives  with  respect  to  translational  temperature  are  stored  as  TYPE 
=  1  while  the  derivatives  with  respect  to  electron  temperature  are  stored  in  the  slot  for  TYPE 
=  2.  NPEEDOT  also  returns  the  derivative  of  tii*  with  respect  to  species  concentrations.  The 
user  should  call  NPEEDOT  even  for  Model  2. 

The  constants  ^"d  in  Eqs.  (58)  to  (60)  are  set  by  NPCHEMIN.  The  effective 
cross  section  a  in  Eq.  (57)  is  computed  in  Subroutine  NPSEET  with  quadratic  fits  taken  from 
Ref.  9. 

The  energy  lost  by  the  free  electron  in  exciting  vibration  is  calculated  in  Subroutine  NPSVE 
(Table  10).  As  discussed  in  Sec.  2.5.2,  NPEEDOT  will  not  call  NPSVE  when  the  vibrational 
and  electron-electronic  energies  are  pooled  (Model  2). 

The  stoichiometric  product  in  Eq.  (61)  is  evaluated  in  NPCHEMIN  and  a  logical  variable, 
IMPACT,  is  set  to  .TRUE,  if  the  product  equals  1. 

3.6  THERMAL  PROPERTIES 

Dimensional  values  {J/kg-mole  or  J/kg-mole/K)  for  the  enthalpy  and  constant  pressure 
molar  heat  capacity  of  the  individual  species  are  obtained  by  calling  Subroutine  NPTHERM 
using  the  calling  sequence  of  Table  IS. 

For  the  one-temperature  model  (Model  1),  NPTABLE  calculates  the  dimensionless  heat 
capacity  and  dimensionless  enthalpy  for  each  species  using  the  polynomial  curve  fits  of 
Eqs.  (62)-(64). 
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Table  14.  Calling  Sequence  for  Subroutine  NPEEDOT 

HPEEDOT (rr , TB , GAN , E vs , CVSV . LOU , IP . I VC » DERIV . EEDOT , EEOOTTH , EBOaTC ) 

IVPOT 


TT(IP) 

T 

Trenslatlonal  teiqiar&tnra  at 
point  n 

(K) 

TE(VP) 

Te 

Electron-el ectronle 

temperature 

(K) 

GAM(MXVS,HP) 

7* 

Concentration  ol  species  s 

at  n 

(kg-moles/m^) 

EVS(MXHS,BP} 

Vilirational  energy  of 
species  s  at  n 

(J/kg-nols) 

CVSV(MXIS,IP) 

Vibrational  specific  beat  of 
species  s  at  n 

(J/kg-nole/K) 

LOV.MP.IXC 

— 

Do-loop  controls 

— 

DERIV 

— 

Flag  to  control  calculation 
of  derivatives 

FROM  COMNON 

NT 

— 

lunber  of  temperature  types 

— 

NS 

— 

lumber  of  species 

— 

NV 

— 

lumber  of  nonequilibrium 
vibrators 

— 

IE 

lumber  of  elements 

OUTPUT 

EEDOT(NP) 

Rate  of  change  of  electronic 
energy 

(U/m®) 

EEDOTTN(HXBT,NP) 

duiljdTn, 

Derivative  of  EEDOT  urt. 
teoperature  of  type  m 

(W/m®/K) 

EEOOTGCHXIS.BP) 

dCjlJd'ii 

Derivative  of  EEDQT  vrt. 
concentration  of  species  s 

(V/kg-mole) 
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Table  15.  Calling  Sequence  for  Subroutine  NPTHERM 
SUBROOTIBE  HPTHEIU!(TT , TV ,TB,LOW , IP , IHC , NODE , IS ,HSV , HSE, CPS , CPSV , CPSE) 


TT(IP) 

TVCMXIV.IP) 

TECIP) 

LOW.IP.IBC 

NODE 

NS 

HV 

HS(HXNS,NP) 

ESV(MXBS.VP) 

NSE(MXKS.NP) 

CPS<HXVS.NP) 

CPSVCNXHS,RP) 

CPS£(MXNS.NP) 


T 

INPUT 

Translational  tenperatnre  at 

(K) 

point  n 

Vibrational  tanperature  ol 

(K) 

T, 

species  s 

Electron- electronic 

(K) 

— 

temperature  at  n 

Do-loop  controls 

_ 

1,2, 3, 4 

Thermal  model  delinition 

— 

(Sec.  2.2) 

FROM  COMMON 

Nunber  ol  species 

— 

Number  ol  vibrators 

— 

h. 

OUTPUT 

Total  species  enthalpy 

(JAg-mole) 

K 

(M0DEL=1) 

Translation-rot  at ion 
part(HODEL>l) 

Vibrational  part  ol  species 

(J/kg-mole) 

hi 

energy 

Includes  electron 

contribution  (K0DEL=2) 
Electronic  part  ol  species 

(J/kg-uole) 

Cp,s 

enthalpy 

Species  molar  heat  capacity 

(JAg-mole/K) 

(MODEL-l) 

Translation-rotation  part 
(MaDEL>l) 

VibrationiLl  part  ol  specilic 

(JAg-mole/K) 

Cl. 

heat 

Electronic  part  ol  specilic 

(JAg-mole/X) 

heat 

MOTES 

1.  Kinetic  energy  el  the  electron  is  listed  as  a  vibration-electronic  energy  lor 
KGDEL=2  and  as  electronic  excitation  energy  lor  H0DEI.=3  or  MODELS. 

2.  Enthalpy  and  energy  ol  internal  nodes  are  equal. 
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IIPUT 

E(HF) 

c™ 

Kixture  energy  of  type  m 

(j/m^) 

GAM(MXHS,HP} 

7* 

Concentration  of  species  s 

(kg-moles/m' 

WHICH 

— 

Flag  to  select  eblcb 
temperature  Is  to  be 
calculated 

LOW.HP.IMC 

— 

Do-loop  controls 

— 

T(HP) 

Initial  guess  for  tbe 
teiq>eratnre 

FROM  COMNOI 

(K) 

NS 

— 

Humber  of  species 

— 

HV 

— 

Humber  of  nonequilibrium 
vibrators 

— 

model 

1.2, 3, 4 

Thermal  model  definition 
(See.  2.2) 

lOIIZED 

— 

If  .TRUE.,  mixture  is 

ionized. 

— 

EFO(NXNS) 

Heat  of  formation  of  species 
8  at  O'*  K 

(J/kg-Bole) 

ETIKOT(iaiS) 

lumber  of 

translation-rotation  modes 

OUTPUT 

T(HP) 

Tm 

Converged  value  of 
tenperaturo  of  type  WHICH 

HOTS 

(K) 

For  MODEL=2,  use  WHICH  -2  for  vibrational  electronic  tenqperatnre.  For  HODEL^S,  use  WHICH 
=2  for  electron-electronic  temperature  and  WHICH  *3  for  vibrational  temperature. 
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In  a  multitemperature  thermal  model,  it  is  necessary  to  separate  the  contributions  of  the 
internal  modes  to  the  species  enthalpy  and  heat  capacity.  This  separation,  which  is  performed 
by  Subroutine  NPTHERM,  is  possible  since  translation  and  rotation  (assumed  fully  excited) 
contribute  a  constant  3.571  to  the  heat  capacity.  If  the  two-temperature  model  (Model  2) 
is  selected,  the  vibrational  portions  of  the  enthalpy  and  specific  heat,  and 
for  molecular  species  s  arc  calculated  by  first  computing  and  Cp^  from  the  curve  fits  using 
7’1/and  then  subtracting  3.5717’^ and  3.57^  respectively.  The  heat  of  formation,  is  also 
subtracted  from  Next,  and  Cp  5  are  simply  set  equal  to  3. 5717' and  3.571for  diatomic 
species  and  2.5717’  and  2.571  for  monatomic  species.  In  both  cases,  the  heat  of  formation 
is  included  with  If  the  gas  is  ionized,  NPTHERM  sets  Ajfand  equal  to  l.SKTy  and 
2.5.71.  If  a  separate  electron  temperature  is  assumed  (Model  a  3),  the  procedure  is  to  evaluate 
Aj  and  Cp  from  the  curve  fits  at  and  the  electronic  contribution  is  calculated  using  7'^, 
in  Subroutine  NPELPART.  Next,  A;  and  are  corrected  as  described  above  to  get  the 
proper  vibrational  contributions.  The  translation-rotation  part  and  the  electronic  part  are 
then  evaluated  at  the  appropriate  temperatures  and  returned  in  separate  arrays. 


The  contributions  of  the  bound  electronic  states  to  the  species  heat  capacity  and  enthalpy 
are  (Ref.  36) 


- 2 - ^  V’jexpC-^,)] 


(69) 


,  SLEVt 
^  t=i 


(70) 


where  =  ^(/kT.  The  electronic  partition  function  is 

NLEV, 

z=  YL 

where  N  LEV^  is  the  number  of  electronic  energy  levels  for  species  s,  g,  is  the  degeneracy 
of  the  fth  electronic  level,  and  is  its  energy. 

3.6.1  Temperature  Calculations 


Since  the  energy  and  heat  capacity  are,  in  general,  nonlinear  functions  of  temperature, 
an  iterative  procedure  is  necessary  to  determine  the  temperature  from  the  energy  or  enthalpy. 
Subroutines  NPEFINDT  and  NPHFINDT  are  provided  for  the  calculation  of  temperature 
from  the  mixture  energy  density  E  and  mixture  enthalpy  density  E  in  J/rr?,  respectively. 
The  calling  sequence  for  Subroutine  NPEFINDT  is  shown  in  Table  16.  The  use  of  NPHFINDT 
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is  identical  except  that  the  input  is  enthalpy  per  unit  voiume  of  the  mixture  rather  than  thermal 
energy  density.  Since  Newton’s  method  is  used  for  the  iteration,  convergence  is  generally 
rapid,  especially  when  the  initial  temperature  estimate  is  good.  The  number  of  iterations  is 
currently  limited  to  ten. 

3.6.2  Entropy 

Some  reacting  gas  codes  (especially  those  which  assume  equilibrium  chemistry)  use  entropy 
as  a  variable.  Subroutine  NPENTRPY  can  be  used  to  calculate  the  species  and  the  mixture 
entropy  for  a  given  translational  temperature.  The  calling  sequence  is  shown  in  Table  17. 

Table  17.  Calling  Sequence  for  Subroutine  NPENTRPY 
SUBKOtlTIIE  IPEITRPY(TT.LOW,iP,IBC,CHI,S,SIIIK) 

IIPUT 


TT(»P) 

T 

Translational  tanperatnra 
array 

(K) 

LOtf.KP.IHC 

— 

Do-loop  controls 

— 

CH1(NXVS,VP} 

X- 

Mole  fraction  of  species  s 

OUTPUT 

— 

S(KIHS,IP) 

— 

Dinens ionless  entropy  of 

species  e 

— 

SNIZ(VP) 

7 

Entropy  of  mixture 

(J/kg^ole/K) 

3.7  TRANSPORT  PROPERTIES 

The  coefficients  of  viscosity  which  enter  the  right-hand  side  of  the  momentum  conservation 
equations,  Eq.  (4),  and  the  thermal  conductivities  which  enter  the  energy  conservation 
equations,  Eqs.  (5)  and  (6),  are  calculated  by  calling  Subroutine  NPTRANS  (Table  18).  In 
order  to  keep  NPTRANS  applicable  to  any  thermal  model,  a  generalized  internal  heat  capacity 
CPI  and  a  corresponding  internal  thermal  conductivity  TKI,  as  in  Eq.  (23),  have  been 
introduced.  Subroutine  NPRETRN  acts  as  a  translator  between  the  notation  of  NPTHERM 
and  that  of  NPTRANS.  Since  it  is  a  translator,  it  must  be  called  between  NPTHERM  and 
NPTRANS.  NPRETRN  determines  the  number  of  internal  energy  storage  modes,  NI,  from 
the  thermal  model  and  the  number  of  nonequilibrium  vibrators  N^,  and  copies  the  heat 
capacities  output  by  NPTHERM  into  the  array  CPI.  Subroutine  NPRETRN  as  written  assumes 
that  the  temperature  of  the  free  electron  is  equal  to  that  of  the  electronic  states.  As  the  number 


53 


AEDC-TR-93-20 


and  sophistication  of  thermal  models  increase,  this  subroutine  will  need  to  be  rewritten.  Usage 
of  NPRETRN  is  detailed  in  Table  19. 

If  the  gas  is  not  ionized,  the  effective  diffusion  coefficients  which  are  returned  by 
NPTRANS  can  be  used  in  Pick’s  law,  Eq.  (15),  to  calculate  the  species  diffusion  velocities. 
If  there  is  only  one  heavy  ion,  the  user  should  employ  the  species  mobilities  returned  by 
NPTRANS  to  calculate  the  ambipolar  diffusion  coefficient  defined  by  Eq.  (21)  for  use 
as  the  effective  diffusion  coefficient  in  Pick’s  law.  If  there  is  more  than  one  heavy  ion  in 
the  chemical  model,  Subroutine  NPOLAR  (Table  20)  must  be  called  to  evaluate  the  species 
diffusion  velocities. 


3.7.1  Transport  Properties  of  Pure  Species 


As  written,  Eq.  (22)  is  applicable  to  a  single-temperature  thermal  model  (Model  1).  For 
multi-temperature  thermal  models  (Model  >  1),  a  separate  coefficient 


M, 


defined  in  Eq.  (23)  is  calculated  for  each  nonequilibrium  mode. 


Note  that  the  transport  properties  of  the  heavy  species  are  evaluated  with  the  translational 
temperature  and  that  rfg  and  are  evaluated  with  the  translational  temperature  of  the  free 
electron  Tg.  For  Model  I,  the  user  must  set  Tg  =  T. 

The  curve-fit  coefficients  and  d/^  in  Eqs.  (65)  and  (66)  for  the  eleven  air  species  have 
been  adapted  from  Ref.  42.  Low-temperature  curve-fit  coefficients  for  species  other  than 
the  components  of  heated  air  have  been  obtained  from  a  variety  of  sources  (Refs.  19  and  41). 

3.7.2  Mixture  Transport  Properties 


The  overall  viscosity  of  the  mixture  is  calculated  by  the  semi-empirical  method  of  Armaly 
and  Sutton  (Ref.  22).  A  similar  approximation,  also  due  to  Armaly  and  Sutton  (Ref.  23), 
is  used  to  obtain  the  individual  parts  of  the  mixture  thermal  conductivity,  for  Model 
<  4.  The  methods  proposed  by  Armaly  and  Sutton  were  adopted  in  NEQPAK  because  when 
NEQPAK  was  under  development  they  were  the  only  models  available  that  specifically  took 
ionization  into  account.  When  varies  from  species  to  species  (Model  4),  it  is 

impractical  to  define  a  mixture  conductivity  for  that  particular  internal  energy  mode.  Hence, 
for  Model  4,  the  contributions  of  species  vibrational  modes  are  returned  as  separate  entries 
in  the  TKI  array. 
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Table  18.  Calling  Sequence  for  Subroutine  NPTRANS 
SUBROUTIIIE  IPTRAISCLOV .IP , IlC . TT ,TB , CPT .CPI .GAN ,D .VIS .NOB . TKT ,TKI , TKE) 


IBPUT 


LOW.IP.IBC 

— 

Do-loop  controls 

TT(IP) 

T 

Heavy  particle  translational 
tenperatnra 

(X> 

TE(VP) 

Te 

Electron  truslational 
tenperatnra 

(K) 

CPT(NXIS,IP) 

Species  translational  heat 
capacity  or  total  heat 
capacity  (NODEL=  1) 

(J/kg-nole/K) 

CPKHXHI.HXIS.IP) 

Species  internal  heat 
capacity 

(J/kg-nole/K) 

GAK(MXHS.IP) 

It 

Concentration  of  species  a 

FROM  COHHOI 

(kg-Boles/m^} 

MI 

— 

lunber  of  nonequilibrinn 
internal  energy  modes 

— 

MODEL 

1.2.3.4 

Themal  node!  definition 
(Sec.  2.2> 

OUTPUT 

DCNXIS.XP) 

D, 

Effective  species 
diffusivity 

<•’/■) 

NOB(HXMS.VP) 

fit 

Ion  mobility 

(m^/s/volt) 

VIS(IP) 

V 

Mixture  viscosity 

(I-s/n*) 

TKT(MP) 

A 

Themal  conductivity  of 
mixture  (H00EL«1)  or  frozen 
aixture  conductivity  (MODEL 
>  1) 

(W/nA) 

TKI(NXII.NP) 

A' 

Internal  conductivity  of 
aixture 

<W/n/K) 

TKEdP) 

Ae 

Themal  conductivity  of 
electrons 

(V/m/K) 
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Table  19.  Calling  Sequence  for  Subroutine  NPRETRN 

SUBROUTIIE  HPRETRKCPS.CPSV ,CPSE .LOW , IP , XIC , CPT, CPI ) 

riPUT 


CPS(NXIS.IP) 

Cv,. 

Species  beat  capacity 

(JAg-mole/K) 

CPSV(NXIS,VP) 

Vibrational  heat  capacity  ot 
species  a 

(J/kg-nole/K) 

CPSE(HXKS,IP) 

Species  electronic  heat 
capacity 

(J/kg-mole/X) 

LOW, IP, lie 

— 

Do-loop  controls 

— 

OUTPUT 

CPT(KXWS,IP) 

Cp.. 

Translational  or  total  heat 
capacity  of  species  s 

(J/kg-noleA) 

CPI(HXVI,HXIS,IP) 

IntameJ.  heat  capacity  of 
node  m  of  species  s 

(JAg-nole/K) 

TO  COMHOI 

— 

lumber  of  noneqnllibrium 
internal  energy  modes 

— 
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Table  20.  Calling  Sequence  for  Subroutine  NPOLAR 
SUBROUTIIE  IP0LAR{G11*,D,MU,G,L0W,W,IIC.V) 

ISPUT 


GAH(KXirS,IP) 

7a 

Concentration  ot  species  s 

(kg-moles/m^) 

D(MXHS,irP) 

D, 

Effective  diflasion 
coefficient  of  species  s 

(«V«) 

MU(HXHS,NP) 

Mobility  of  species  s 

(m^/volt/s) 

G(NXVS.n>) 

VXa 

Gradient  of  mole-fraction  of 
species  s 

(M“») 

LOW.XP.IXC 

Do-loop  controls 

FROM  CDNMOl 

iLPZJ(MXVE,mS) 

“m 

Mumber  of  atoms  of  kind  "i" 
in  species  "j" 

OUTPUT 

V(MXXS,HP) 

«a 

Diffusion  velocity  of 

species  s 

(kg-moles/m^) 

57 


AeDC-TR-93-20 


3.8  UTILITY  ROUTINES 

NEQPAK  provides  two  utility  routines,  NPFRCON  and  NPTOCON  (Tables  21  and  22) 
to  convert  composition  units  to  and  from  the  kg-moles/m^  used  in  the  computational 
modules.  Subroutine  NPWEIGHT  (Table  23)  calculates  the  molecular  weight  of  the  mixture. 

Table  21.  Calling  Sequence  for  Subroutine  NPFRCON 

SUBROUTllB  IPFllCOV(ARiUY,LOV,IP,INC,TO.DEMS,UT.OUT) 

Thio  subroutine  converts  composition  or  rate  units  from  kg-moLes/n^  to  the  units 
specilied  by  the  character  variable  TO. 

INPUT 


ARRAY(KX1IS,NP) 

— 

Array  of  compos it ions  or 

— 

rates 

LOV.HP.XIC 

— 

Do-loop  controls 

— 

TO 

_ 

Units  of  "OUT".  Use 

— 

*aole/moIa' , ’gram/gram* , 

‘mole/gran* 

or  'gram/cn.m'.  Enclose 
single  qnotesC') 

in 

DERS(NP) 

P 

Mixture  mass  density 

(kg/m®) 

WTCNP) 

M 

Mixture  molecular  weight 

(kg/kg-nole) 

OUTPUT 

□UT(BP) 

— 

Output  array 

— 
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Table  22.  Calling  Sequence  for  Subroutine  NPTOCON 

SUBKOUTIVB  HPTOCOK ARUT .  LOU ,  VP ,  IKC ,  FROM .  DEVS .  VT .  OUT) 

This  subroutine  converts  composition  or  rate  units  from  the  units  specified  by  the 
character  variable  FROM  to  hg-moles/m^ . 


IVPUT 


ARRAY(HXVS.VP) 

— 

Array  of  compositions  or 

rates 

LOV.VF,IVC 

— 

Do-loop  controls 

FRDK 

ttolts  of  "ARRAY".  Use  — 

*nols/aole ’ , *gram/gram* , 'mole/gram' , 
*grBm/cu.m’ , *mole-fraction* 
or  'mass-fraction*.  Use 

lover  case  and  enclose  in 
single  quotesC*) 

DEVS(VP) 

P 

Hizture  mass  density 

(Rg/m^) 

UT(HP) 

M 

Mixture  molecular  weight 

FROM  CQHHOR 

(kg/kg-mole) 

EHI(KXMS) 

Ms 

Species  molecular  weight 

OUTPUT 

(kg/kg-mole) 

auT(vp) 

Output  array 

(kg-mole/m^ ) 
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Table  23.  Calling  Sequence  for  Subroutine  NPWEIGHT 

SUBaaUTIHE  IPVEIGHT(C , LOW . HP . IHC .UVITS . WT) 

IHPUT 

Array  ol  compoBltions  — 

Do-loop  controls  — 

Units  oi  ''C.  Use  — 

'boIo/doIs' , ’gran/gran* , 'mole/gran' , 
'gran/cn.n' .  Use  lower  case 
and  enclose  in  single 
qnotesC'} 

OUTPUT 

Mixture  molecular  weight  (kg/kg-mole) 

3.9  EQUILIBRIUM  COMPOSITION 


CCMXITS.BP} 

LOW, HP. lie  — 

UIITS  — 

WT(HP)  M 


Interpolation  for  a  particular  temperature  and  pressure  within  the  table  of  equilibrium 
compositions  prepared  by  Subroutine  NPREPEQ  (Table  6)  is  performed  by  Subroutine 
NPEQCOMP  (Table  24).  Note  that  Fortran  unit  lEQ  must  have  been  opened  as  a  direct 
access  device  before  calling  NPEQCOMP  for  the  first  time.  The  record  length  is  8*(MXNS 
f  4). 


I 


I 

I 

I 

J 


60 


AEDC-TR-93-20 


Table  24.  Calling  Sequence  for  Subroutine  NPEQCOMP 
SUBROUTZirE  IPBqC0HP(TVALL,PVALI..IEq,O(rr) 


INPUT 

TWALL 

— 

Temperature 

(K> 

PUALL 

— 

Pressure 

(Pa) 

lEQ 

— 

Fortran  unit  number  for 
input 

— 

FROM  COKHOir 

HEQT  —  lumber  ot  aq;ailibriiim  — 

tenperatureB 

NEQP  —  Vumber  ot  equilibrium  — 

preasurea 

DELT  —  Teiqterature  iucrement  (K) 

DELP  —  Pressure  Increment  (Pa) 


OUTPUT 

OOT(MXI»S+2)  —  Output  array  — 

VOTE 

The  Arst  element  of  the  OUT  array  contains  the  density  of  the  mixture  in  kgfm^^  while  the  second 
clement  contains  the  mixture  frozen  sound  speed  in  in  m/s.  The  remaining  N,  positions  axe  iilled 
with  the  species  concentrations  in  kg  —  moles /m^.  The  compositions  axe  listed  in  the  order  in 
which  they  are  returned  in  the  FSCQNP  array  by  Subroutine  NPCHEMIN  (Sec.  3.2.1). 
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4.0  CONCLUDING  REMARKS 

NEQPAK  is  a  set  of  Fortran  77  routines  and  databases  which  are  useful  in  the  development 
and  application  of  both  reacting  flow  solvers  and  static  kinetics  programs.  It  is  unique  in 
its  inclusion  of  thermal  nonequilibrium  effects  and  in  the  range  of  its  databases. 

■  The  thermal  properties  database  covers  the  range  from  100  to  35,000  K. 

•  The  transport  properties  database  extends  from  300  to  30,000  K. 

•  Ambipolar  diffusion  effects  may  be  evaluated. 

•  The  effects  of  thermal  nonequilibrium  on  the  thermal  conductivity  may  also 
be  evaluated. 

•  A  utility  for  correlating  user-supplied  reaction  data  with  the  material  properties 
databases  is  provided. 

•  The  interactions  between  vibrational  and  chemical  relaxation  may  be  investigated 
by  using  multiple-temperature  reaction  models  or  by  application  of  coupling 
factors  suggested  by  contemporary  theories. 

•  A  module  for  calculating  the  equilibrium  chemical  composition  as  a  function 
of  pressure  and  temperature  is  also  included.  This  module  may  be  of  use  in 
initializing  flow  calculations  or  in  evaluating  solid  surface  boundary  conditions. 

The  modularity  and  flexibility  of  NEQPAK  invite  extension  of  its  capabilities  in  many 
directions.  The  most  pressing  and  practicable  enhancements  are 

•  extension  of  the  transport  database  to  the  very  low  temperatures  encountered 
in  wind  tunnel  testing, 

•  the  development  of  a  surface  chemistry  capability,  and 

•  the  capability  to  calculate  the  energy  lost  through  equilibrium  radiation. 
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APPENDIX  A 

SPECIES  EQUATION  WITH  AMBIPOLAR  DIFFUSION 

This  appendix  derives  the  form,  Eq.  (18),  of  the  species  conservation  equation  which  applies 
to  a  plasma  containing  several  positively  and  negatively  charged  ionic  species.  In  the  following 
derivation,  F  =  7/  and  all  summations  have  the  limits  1  and  N,. 

Substitution  of  Eq.  (7.3-27)  into  Eq.  (7.4-3)  of  Ref.  7  and  neglect  of  the  thermal  and 
pressure  diffusion  terms  yields  the  following  equation  for  the  diffusive  mass  flux,  which  appears 
on  the  right-hand  side  of  Eq.  (3). 

i=i  P  ^  k  (A-1) 

The  units  in  Eq.  (A-1)  are  consistent  (kg-mole-m^/s)  if  ^  is  a  force  per  kg-mole  and  the  units 
of  the  multicomponent  diffusion  coefficient  are  m^/s.  Note  that  D,,  =  T>^  only  if 
=  2.  Let 


=  -ae,,^feS  (A-2) 

where  e  is  the  elementary  charge,  E  is  the  electric  field  strength  (N/C)  and  is  Avogadro’s 
number.  The  number  of  electrons  per  particle  of  species  s  is  Note  that  a  for  a  positive 
ion  is  negative  since  one  or  more  electrons  have  been  removed.  Substitution  of  Eq.  (A-2) 
into  the  flux  equation,  Eq.  (A-1),  yields 


(A-3) 


From  Eq.  (A-3)  it  is  clear  that  the  diffusive  flux  contains  a  term,  known  as  the  “drift  velocity,” 
which  is  proportional  to  the  electric  field  strength!?.  The  constant  of  proportionality  is  known 
as  the  “ion  mobility,”  For  a  neutral  plasma,  S*  =  0  and 

7  ?  +  -yn/ja^jS] 

It  was  shown  in  Table  1  of  Sec.  2.3  that  in  the  absence  of  ionization,  Fick’s  law 


=  -TD^Vx,  (A-5) 

is  a  good  approximation  to  Eq.  (A-1)  if  the  effective  diffusion  coefficient  is  appropriately 
defined.  Now,  if  it  is  assumed  that  the  right-hand  side  of  Eq.  (A-5)  continues  to  be  a  reasonable 
approximation  to  the  first  term  of  Eq.  (A-4),  even  in  the  presence  of  an  electric  field,  and 
an  ion  mobility  is  defined  as 
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_  T^eAf^j 

”  PP  7*tte,.  (A-6) 

then  the  species  conservation  equations  become 


^  +  V  .  (7.Cr)  =  «).  +  V  •  [rP.Vx.  - 


(A-7) 


Equation  (A-7)  is  still  general.  To  specialize  Eq.  (A-7)  for  ambipolar  diffusion,  multiply  both 
sides  by  a,  ^  and  sum  over  The  result  is 


dt 


53  +  ^  ■  [r  ]C 


(A-8) 


Conservation  of  charge  requires  that 

N» 


{A-9) 


Since  the  plasma  is  neutral, 

JVj 

53  =  0 

*=1  (A-10) 

Then  Eq.  (A-8)  becomes 


V  •  [r  5]  a^^y.DtVxs  -  $3  “2, */*»7»^l  =  0 

»  •  (A-11) 

In  general,  the  electric  field  has  both  applied  and  induced  components,  E  =  +  Ej.  If 

there  is  no  applied  field,  i.e.,  E  is  entirely  due  to  the  equilibrium  between  the  electrostatic 
and  collision  forces,  then 


&i  =  -Y\ct,,,D,Vxs 

where 


(A-12) 


o  =  (A-13) 

J 

WhenE  =  Sj,  defined  by  Eq.  (A-12),  is  substituted  into  Eq.  (A-7),  the  ambipolar  species 
equation,  Eq.  (18),  is  recovered,  namely 


^  +  V .  (7,t?)  =  li.  +  V .  [TCe.Vxs  -  2:^  E  o.jH,  Vx,)) 


(A-14) 
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Kinetic  theory  attempts  to  calculate  the  effect  of  intermolecular  collisions  on  the  distribution 
of  molecular  population.  The  formulae  usually  quoted  for  thermal  conductivity  apply  only 
to  monatomic  particles  and  must  somehow  be  corrected  for  the  effects  of  internal  energy. 
While  there  are  formal  treatments  for  polyatomic  molecules,  the  formulae  are  cumbersome 
and  are  also  restricted  to  thermal  equilibrium  and  by  various  assumptions  about  relaxation 
rates.  Therefore,  most  authors  simply  apply  some  version  of  the  Eucken  correction.  The 
following  analysis  follows  that  in  Hirschfelder,  Curtiss,  and  Bird  (Ref.  7,  pp.  498-SOl)  with 
some  slight  changes  in  notation.  A  chemically  pure  polyatomic  gas  with  a  set  of  internal 
quantum  numbers  {1}  is  assumed.  The  essential  step  in  the  following  treatment  is  the 
assumption  that  each  possible  quantum  state  of  the  gas  constitutes  a  separate  species.  Thus, 
the  single  polyatomic  gas  is  replaced,  in  thought,  with  a  mixture  of  monatomic  species  to 
which  the  results  of  ordinary  kinetic  theory  apply. 

One  of  the  general  results  of  the  Chapman-Cowling-Enskog  (Ref.  7)  solution  of  the 
Boltzmann  equation  is  that  the  flux  of  energy  across  some  reference  surface  is 

f  =  -A°VT  +  2  »i(5*r/2  +  Ci)«;  (B-1) 

t 

where  e-,  is  the  energy  per  molecule  of  the  quantum  state  ;  and 

A®  =  5f^CV/2  =  (15*/4m)»/ 

is  the  thermal  conductivity  due  to  translation.  The  viscosity  is  1;,  the  mass  of  the  molecule 
is  m,  and  is  the  number  density  of  the  molecules  which  are  in  the  fth  quantum  state.  The 
quantities  n  =  and  X,  =  n/n  are  also  used. 

It  is  a  good  approximation  that  the  diffusion  velocity  is  independent  of  the  internal  quantum 
state.  Then,  since  the  masses  of  all  these  species  are  the  same,  the  diffusion  velocity  of  species 
/  is 

=  -XT^V^Xi  (8-2) 

where  "D  is  the  self-diffusion  coefficient.  It  is  possible  to  rewrite  Eq.  (B-1)  as 

?  =  -  A®Vr  +  2.bkT  ^ -H  (B-3) 

t 

Since  the  masses  of  all  the  excited  species  are  equal,  the  deBnition  of  diffusion  velocity  implies 
that  Ln/n-Ui  =  0  =  in  this  case.  Substitution  of  Eq.  (B-2)  into  Eq.  (B-3)  yields 
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(B-4) 


Equation  (B-4)  is  valid  for  both  thermal  equilibrium  and  thermal  nonequilibrium.  For  thermal 
equilibrium,  the  distribution  of  population  among  the  internal  states  is  Maxwellian  with  T 
as  the  temperature  parameter.  Since  the  temperature  is  now  defined,  application  of  the  chain 
rule  for  derivatives  and  recognition  that  e,  is  a  function  only  of  /  yield 


?'*=  -AoVr  -  =  -AoVT  - 


(B-5) 


in  which  Of,=  e,  dX/dT  is  the  heat  capacity  of  the  rth  species.  Equation  (B-5)  is  obtained 
by  defining  Cp  a  Cy+  H  and  recognizing  that  the  Cp  returned  by  Subroutine  NPTHERM 
contains  a  contribution  2RT/2  from  the  translational  motion  which  is  included  in  the  Xg  term 
and  a  contribution  KT  from  the  pV-viorV.  part  of  the  enthalpy. 


In  a  multitemperature  model,  it  is  assumed  that  the  first  quantum  states  correspond 
to  rotation  with  a  Boltzmann  temperature  that  states  1,.  <  i  £  correspond  to  vibration 
and  are  distributed  according  to  r„,  while  the  remaining  states  <  i  £  Ig  correspond  to 
electronic  excitation  with  a  Boltzmann  temperature  Tg.  The  energy  flux  can  then  be  written 


ir 


g=  -AoVT  -  nD^C\,VTr  +  •••]=  -AoVT  - 


i=l 


where  V  =  nDC*^&  the  thermal  conductivity  of  the  /th  internal  mode. 
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Variable  Type 

ALPIJ(MXNE,MXNS)  R 

AK(MXNR),..-,CK(MXNR)  R 

DK(MXNR),...,FK(MXNR)  R 


GK(MXNR),...,UK(MXNR)  R 
B(MXNS,MXNS)  R 

CHI1,CHI2,CHI3  R 

CONTK(MXNS,MXNS)  R 

CONV(MXNS,MXNS)  R 

COUPLE  L 

CTYPE(MXNR)  I 

C  ID,. . .  ,C5D(MXNS,MXNS)  R 
DEGEN(7,MXNS)  R 

DISOC(MXNS)  R 

EFF(MXNS,MXNR)  R 

ELEV(7,MXNS)  R 

EMI(MXNS)  R 

EMIQ(MXNS)  R 

ETNROT(MXNS)  R 

EVCOEF(5,MXNS)  R 

HFO(MXNS)  R 

IMPACT(MXNR)  I 

lON(MXNS)  R 

IONIZED  L 

LINK(MXNS)  I 

LION  I 

METHOD  L 

MODEL  I 

NC  I 

NE  I 

NI  I 

NR  I 

NS  I 

NT  I 

NTBR(MXNS)  I 

NTHD  I 

NU(MXNR)  R 

NUP(MXNR,MXNS)  I 

NUR(MXNR,MXNS)  I 


Definition 

Number  of  atoms  of  element  i  in  species  j 
Parameters  in  forward  reaction  rate 
A,T®-exp(-Cr/T),Eq.(28) 

Parameters  in  reverse  reaction  rate 
Drr®-exp(-n/T),  Eq.  (67) 

Parameters  in  Sec.  3.3.5 

Constant  set  in  NPVISCIN.  Used  in  NPSUTTON 

Parameters  fi , .  -  • ,  fs  in  Sec.  3.5 

Constant  defined  in  NPVKCIN.Used  in  NPSUTTON 

Constant  defined  in  NPVISCIN.  Used  in  NPARMALY 

Switch  used  in  NPSOURCE  to  caU  NPCOUPLE 

Coupling  type  for  reaction  r,  Sec.3.3.5 

Coefficients  of  polynomial  fit  to  binary  diffusivity 

Degeneracy  of  electronic  energy  level 

Dissociation  energy 

Catalytic  efficiency 

Energy  of  electronic  level 

Species  molecular  weight 

Fourth  root  of  spedes  molecular  weight 

Dimensionless  enthalpy  of  rigid  rotator 

Fit  coefficients  used  in  NPSVE 

Heat  of  formation 

Marker  for  impact  ionization  reaction 
Ionization  energy  of  spedes 
Flags  an  ionized  mixture 

Links  virtual  molecule  to  real  counterpart.  Used  in  NPCOUPLE. 

Number  of  heavy  ions 

Used  in  NPSOURCE  to  call  NPGIBBS 

Distinguishes  thermodynamic  models 

Maximum  value  of  CTYPE(IR) 

Number  of  elements 

Number  of  internal  energy  modes.  Used  in  NPTRANS 

Number  of  reactions 

Number  of  spedes 

Number  of  temperature  types 

Number  of  thermo- fit  polynomials 

Number  of  three- body  reactions 

Change  in  number  of  moles  in  reaction  IR 

Stoichiometric  coefficient  of  product 

Stoichiometric  coefficient  of  reactant 
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REV(MXNR) 

SIGMA(3,MXNS) 

SLJ(MXNS,MXNS) 

SUA(MXNS) 

SUBB(MXNS) 

TAB(7,MXBR,MXNS) 

TAUC0N(4,MXNS,MXNS) 

TBR(MXNS,MXBR) 

THET(MXNS) 

TTYPEF(MXNR) 

TTYPER(MXNR) 

ZER0(MXNS,MXNS) 

ZROT(MXNS) 


I  flags  method  for  computing  reverse  reaction  rate 
R  Fit  coefficients  used  in  NPSEET 
R  Lenard- Jones  collision  diameter 
R  Constant  in  Sutherlsmd  viscosity  formula 
R  Constant  in  Sutherland  viscosity  formula 
R  Thermo-lit  coefficient 
R  Coefficients  of  relaxation  time  fit 
R  Break  point  between  thermo-fit  segments 
R  Characteristic  temperature  for  vibration 
I  Temperature  type  in  forward  direction,  Sec.  3.3.1 
I  Temperature  type  for  reverse  direction 
R  Is  zero  on  diagonal  and  one  elsewhere 
R  rotational  collision  number 
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All  of  the  NEQPAK  modules  are  listed  in  Table  D-1  below.  Those  subroutines  which 
must  be  called  by  the  user  are  marked  by  the  word  “FLOW”  in  the  “Called  by”  column. 
The  calling  sequence  of  these  modules  is  presented  in  the  text  in  the  indicated  tables. 


Table  D-1.  NEQPAK  Modules 


Module 

Called  by 

Function 

Table 

gamin 

gsei^cf 

helps  calculate  incomplete  gamma  function 

gcf 

npsgam 

helps  calculate  incomplete  gamma  function 

gser 

npsgam 

helps  calculate  incomplete  gamma  function 

Itrim 

npchemin 

length  of  character  string  with 

npdifhn 

trailing  blanks  removed 

nparmaJy 

nptrans 

returns  mixture  viscosity 

npait 

npcouple 

returns  partition  function  of  harmonic 

npsvc 

oscillator  truncated  at  energy  ’’a” 

nparta 

npcouple 

derivative  of  npart  with 

npsvc 

respect  to  truncation  energy 

npartt 

npcouple 

derivative  of  npart  with  respect  to  temperature 

npchemin 

FLOW 

chemistry  input 

5 

npcopy 

npchemin 

copies  species  data 

np  couple 

npsource 

returns  vibrational  coupling  factor 

npcphs 

npeqlbrm 

returns  specific  heat  and  enthalpy 

npdifHii 

npchemin 

selects  diffusion  data 

npeedot 

FLOW 

electron  energy  source  term 

14 

npehndt 

FLOW 

finds  temperature  for  a  given  energy 

16 

npdpart 

nptherm 

energy  in  electronic  states 

npentrpy 

FLOW 

returns  species  entropy 

17 

npeqcomp 

FLOW 

interpolates  in  composition  table 

24 

npeqlbrm 

npiepeq 

returns  equilibrium  composition  for  one  p,T 

npevdot 

FLOW 

vibrational  energy  source  term 

9 

npfrcon 

FLOW 

converts  composition  variables 

21 

npgauss 

npeqlbrm 

inverts  matrix 
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Table  D-l.  NEQPAK  Modules  (Concluded) 


Module 


Called  by  |  Function 


I  l^ble 


npgibbs 

npbflndt 

Qplaxin 

npmatrix 

npnopref 

upolar 

nprelax 

nprepeq 

npretrn 

npsec 

npseet 

npseten 

npsgam 

npsouice 

npsutton 

npsvc 

npsve 

npsvv 

nptable 

nptherm 

nptocon 

nptrans 

npvisdn 

npweight 


npsouice 

PLOW 

npchemln 

npeqlbrm 

FLOW 

npevdot 

FLOW 

FLOW 

FLOW 

npeedot 

npeedot 

npiepeq 

npcouple 

FLOW 

nptrans 

FLOW 

FLOW 

npeedot 

FLOW 

nptherm 

npefindt 

FLOW 

FLOW 

FLOW 

npchemln 

FLOW 


returns  Gibbs  free  energy 

finds  temperature  for  a  given  enthalpy 

sets  data  for  relaxation  time 

sets  up  composition  matrix 

returns 

ambipolar  ion  diffusion  velocity  20 

returns  vibrational  relaxation  time  13 


driver  for  chemical  equilibrium  calculation 

sets  up  call  to  nptrans 

loss  of  energy  by  electron  impact 

electron  energy  exchange  by  elastic  collision 

initializes  composition  for  npeqlbrm 

returns  incomplete  gamma  function 

species  source  term 

returns  mixture  thermal  conductivity 

effect  of  chemistry  on  vibrational  energy 

excitation  of  vibration  by  electron 


6 

19 


7 

12 

10 


vibration- vibration  exchange 
dimensionless  enthalpy  and  Cp  for  species  s 


11 


enthalpy  and  Cp  for  species 
composition  variables  to  concentration 
returns  transport  properties  of  mixture 
sets  up  viscosity  data 
calculates  mixture  molecular  weight 


15 

22 

18 

23 


76 


AEDC-TR-93-20 

APPENDIX  E 

ENTERING  REACTION  DATA 

All  reaction  data  are  list  directed  input  and  are  read  by  Fortran  unit  INPR.  Data  items 
are  separated  by  commas.  Omitted  items  are  represented  by  commas.  Items  following  a  slash 
“/”  will  not  be  read. 

The  first  record  of  unit  INPR  is  used  to  specify  the  input  units  of  the  rate  coefficients 
and  the  name  of  the  reaction  set,“RNAME”.  Use  “cu.m/kmol”,  “c.cm.mole”,  “c.cm/part”, 
or  “cu.m/part”  in  lower  case  letters  as  shown  to  designate  the  units.  The  default  option 
is  “c.cm/mole”. 

For  each  reaction  there  is  a  principal  record  which  defines  the  reaction  and  Identifies  the 
temperature  functions  to  be  used  in  the  forward  (TYPEF)  and  the  reverse  (TYPER)  directions 
and  also  identifies  the  coupling  model  (KPL).  The  six  dummy  variables  A,  ....  F  describe 
the  temperature  dependence  of  the  forward  and  reverse  rates  as  described  in  Secs.  3.3.2  and 
3.3.3.  The  meaning  of  the  parameters  G,  H,  and  t/,  which  describe  the  effea  of  nonequilibrium 
vibration  on  the  reaction,  depends  on  the  value  of  KPL  as  described  in  Sec.  3.3.5.  The  format 
for  this  principal  record  is 

REACTION,  TYPEF,  TYPER,  KPL,  A,B,C,D,E,F,G,n,U 

where  REACTION \s  defined  as  CHARACTER*63.  TYPEF  and  TYPER  are  integers.  Species 
names  are  limited  to  eight  characters,  must  begin  with  an  alphabetic  character,  may  not  contain 
embedded  blanks,  and  are  delimited  with  spaces.  NEQPAK  will  stop  and  print  an  error  message 
if  a  species  name  cannot  be  found  in  the  thermochemical  database.  Species  names  in  the 
current  database  contain  only  upper  case  letters.  An  irreversible  reaction  is  indicated  by  using 
a  minus  sign  in  place  of  the  “  =  ”  as  a  separator  between  reactants  and  products. 

The  activation  temperature  in  the  argument  of  the  exponent  is  a  positive  number 
corresponding  to  Eq.  (28).  Care  must  be  taken  when  using  published  tabulations  of  reaction 
rates  as  there  is  no  universally  accepted  standard  for  units  or  for  the  sign  of  C^.,  and  some 
authors  are  careless  about  stating  their  conventions. 

For  three-body  reactions,  indicated  by  the  appearance  of  the  character  “M**  among  the 
reactants,  an  auxiliary  record  is  used  to  input  catalytic  efficiencies.  The  format  of  the  auxiliary 
record  is 

{OBJECT{N),VALUE{N),N  =  1,MXNS), 
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where  OBJECT  is  defined  to  be  CHARACTER*8,  VALUE  is  a  real  number,  and  MXNS 
is  the  maximum  number  of  species  (currently  20)  set  by  the  PARAMETER  statement  in  file 
NPDIMEN. 

A  typical  three-body  reaction  might  be  entered 

*02  N  °  0  +  0  '»■  N*,1,1,0,3.61E18,-1.0.59400..3.01E15,-0.6/ 
*D2*,9.,*0’,2S..’N2*,2.0.'NO+*,0.0,'E-*,0.0/ 

Note  that  stoichiometric  coefficients  are  not  entered,  that  the  letter  “M"  represents  any  third 
body  and  causes  the  auxiliary  record  to  be  read,  and  that  the  rate-controlling  temperature 
is  specified  by  the  TYPEP  value  of  one  —  usually  the  translational  temperature.  Furthermore, 
the  reaction  is  reversible  since  the  separator  between  reactants  and  products  is  an  “  =  ”  sign. 
Also,  D  and  E  ^  0.0  so  that  REV(r),  Sec.  3,3.3,  is  set  automatically  =  2.  The  auxiliary 
fields  specify  catalytic  efficiencies  ^  of  25.0  for  oxygen  atoms,  9.0  for  oxygen  molecules 
O2,  and  2.0  for  nitrogen  molecules  Nj.  All  other  species  have  ^=1.0  except  the  ions  NO* 
and  E~  which  do  not  affect  the  reaction. 

Photochemical  reactions  may  be  specified  by  entering  the  character  string  ‘‘HNU’’  as 
either  a  product  or  a  reactant.  The  OBJECT  and  VALUE  entries  on  the  auxiliary  record 
then  represent  an  identifier  for  the  photon  and  its  wavelength,  respectively.  Needless  to  say, 
this  identifier  must  be  listed  as  a  species  in  the  thermal  properties  database. 

As  long  as  the  rate  equation  fits  into  the  Arrhenius  form,  Subroutine  NPSOURCE  will 
return  appropriate  gain  and  loss  terms.  Since  photochemical  reactions  are  often  essentially 
irreversible,  the  user  must  take  great  care  to  ensure  mass  balance  in  setting  up  the  reaction 
model.  NEQPAK  has  no  mechanism  for  checking  this  balance,  nor  does  it  possess  a  procedure 
for  handling  the  radiant  energy. 
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The  first  three  reactLon  sets  shown  below  were  prepared  to  serve  as  examples  in  Sec.  3.3.5. 
They  are  basically  the  same  as  the  Blottner  set,  except  for  specifying  the  indicated  coupling 
and  using  the  electronic  temperature  where  appropriate.  These  reaaion  sets  were  not  necessarily 
used  by  the  indicated  authors.  The  remaining  data  sets  have  been  translated  to  NEQPAK 
notation  and  set  up  as  Model  1  cases.  The  hydrogen-air  example  is  taken  from  Ref.  52. 

F-1.  PURE  AIR  EXAMPLES 


‘c. cm/mole, hassanV  (Ref.  1 1) 

'02  +  H  -  0  +  a  +  H'. 1,1, 1.3. 61E18,-1. 0.59400. ,3. OIEIS. >0.5/ 
'ND+'.0.0.'E-'.0.0.’02' ,9..’0',2S.,*N2'.2.0/ 

'N2  M  B  H  +  N  4  M',l,l,l,1.92E17,-0.5,1.131ES,l.09E16,-0.6/ 
'N'  ,0.0,'l!rO+',0.0,'E-’,0.0,'N2»,2.5/ 

’N2  4N  =  H  +  ir4ir  ',1,1,1,4.15E22,-1.6,1.131E5,2,32E21,-1.S/ 

'NO  4  M  =  H  4  0  4M',1,1,1,3.97E20,-1.6,7.66E4,1,01E20,-1.5/ 

'KD4»,0.0,'E-',0,0,’D*,20.,*M',20.,'N0*,20.y 

'NO  4  G  •  02  4  IfM,l,0,3.18E9,1.0,1.97E4.9.63E11.0.5.3600./ 

*N2  4  0  -  HQ  4  H  », 1. 1, 0. 6. 75E13, 0.0, 3. 7SE4,1.SE13. 0.00,0.0/ 

‘M  4  0  -  1104  4  E-*,l,2.0,9.03E9,0.5,3.24E4,1.8E19,-i.O,0.0/ 
'END  ' 


‘c.cra/mole,candierV  (Ref.  12) 

'02  4  H  -  0  4  G  4  H', 1,1,0, 3. 61E18.-1.0, 59400., 3.01E1S, -0.6/ 
’02’,9.,'0',25.,'N2',2.0/ 

*N2  4  M  «  H2*  4  M',1,1,2,8.26E16,0.0,25289.,0.,0.,0.,5.0,0.36/ 
'N',0.0,/ 

'N2  4  N  =  N2*  4  N  1.1, 2, 6. 028E16,0.,0.,0.,0.,0. ,5.0,0.36/ 

’N2*  4  M  -  N  4  K  4  M',1,1.0,1.S2E21,-1.51,13121.,0.,0.,0.,6.0.0.35/ 
'N’,0.0/ 

'N2+  4  M  -  g  +  IT  +  M',1,1,0.6.76E21,-1.S1.13121.,0.,0.,0..5.0,0.35/ 
'HQ  4  M  =  K  4  0  4M’,1,1,0,3.97E20,-1.S,7.S6E4.1.01E20.-1.6/ 
'0’,20.,'N’.20.,'N0',20./ 

'NO  4  0  -  02  4  N’, 1, 1. 0, 3. 18E9, 1.0, 1.97E4, 9. 63E11, 0.6,3600./ 

'H2  4  D  -  NO  4  N  M, 1,0,6. 76E13,0.0, 3. 76E4,1.5E13, 0.00,0.0/ 

'N  4  0  »  N04  4  E-M,2,0,9.03E9,0.6,3.24E4.1.8E19,-1.0.0.0/ 

'END  ' 
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‘c.cm/mole.fruhaufV  (Ref.  10) 

'02  +  H  =  0  +  0  +  MM,l,4,3.61E18.-1.0,69400..3.01el5.-O.B,,0. 
'W0+»,0.0,'E-’,0.0,>02'.9..'G',2S..'N2',2.0/ 

'>12  +  H  =  N  +  V  +  MM.1.4,1.92E17,-0.5.1.131E6,1.09el6,-0.5..0 
'N'.O.O.'NQ+'.O.O.'E-' ,0.0,*H2',2.5/ 

'N2-i-lfeN  +  N4^H  M,1.0,4.1BE22.-1.5,1.13iES,2.32E21,-1.5/ 
'NO  +  M  =  M  +  0  +MM,1,0.3.97E20,-1.6,7.56E4,1.01E20,-1.6/ 
’lfD4',0.0.'E-',0.0,'0'.20..'Jf>,20..'HD',20./ 

'NO  +  D  -  02  +  NM,1.3,3.18E9.1.0,1.97E4.9.636ll,0.6.3600.,0.7 
'112  +  0  -  KO  +  N'. 1. 1, 3. 6. 76E13, 0.0, 3. 75E4. 1.5813,0.0, 0.7, 9400, 
'JT  +  0  =  SO*  +  E-M, 2,5,9. 03E9, 0.6, 3. 24E4. 1.8819,-1.0,0.0, 0.7, 
'END  ' 


‘c.cm/raole.blottner’/  (Ref.  51) 

'02  +  H  -  0  +  0  +  M',1,1.0,3.61E18,-1.0,6940.0,,3.01E1S,-0.6/ 
'NO*' ,0.0, 'E-', 0.0/ 

'N2  +  M  =  N  ♦  N  +  H',l,l,0,1.92E17,-0.6,1.131ES,.1.09E16,-0.5,/ 
'N',0.0,'NO+’,0.0,'E-»,O.Q  / 

'N24N-N4K4N  ',l,l,0,4.16E22,-l.B.1.131E6,,2.32E21.-l.S/ 
'BO  +  H  »  H  +  0  +M',1,1,0,3.97E20,-1.6,7.56E4,,1.01E20,-1.6/ 
’B0+',0.0,'E-»,0.0/ 

'MO  +  D  -  02  +  N’,1.1,0,3.18E9,1.0,1.97E4..9.63E11,0.5,3600/ 

'M2  +  0  ■  MO  +  N  M,1,0,6.7513,0.0,3.7SE4,,1.6E13/ 

'M  +  0  -  MG+  +  E-',l,l,0,9.03E9,0.6,.3.24E4,,1.8E19,-1.0/ 

'END'/ 


7.14875.0/ 

.7,28275./ 


,18875./ 

/ 

8.183/ 
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‘c.cin/mole,guptaV  (Ref.  42) 

•02  +  M  »  D  +  0  +  MM,1.0,3.61E18.-1.0,5940.0,,3.01E16.-0.5/ 

•N0+' ,0.0,*E-',0.0/ 

•N2  M  =  N  +  N  +  M>,l,l,0,l.92E17.-0.S,1.131E5,,1.09E16,-0.5,/ 
'N’,0.0,'MO+'»0.0,^E-',0.0  / 

•N2  +  N«K+N  +  N  M,l,0.4.15E22,“1.5,1.131ES,.2.32E21,“i.6/ 

•NO  +  M  »  H  +  a  +M  ’.1.1,0,3.97E20,-1.S,7.56E4..1.01E20.-1.S/ 
•N0+',0.0,’E-\0.0/ 

•NO  +  0  s  02  +  N  '3l.l,0.3.18E9.1.0,1.97E4..9.63Eli,0.5,3600/ 

'N2  +  0  -  NO  +  N  ’,1.1.0.6.7513,0.0,3.75E4.,1.5E13/ 

•N  +  0  -  ND+  +  E-  »,1,1,0,9.03E9.0.5,,3.24E4.,1.8E19.-1.0/ 

•0  +  E-  »  0+  +  E-  +  E-  ’,1.1.0.3.6E31,-2.91.1.58ES,,2.2E40,-4.5  / 

•N  +  E-  =  N+  +  E-  +  E-  M,l,0,l.lE32,-3.14,1.69E5,,2.2E40,-4.5/ 

•0  +  0  *  02+  +  E-  M,1,0,1.6E17.-.98,8.08E4.,8.02E21.-1.5/ 

'0  +  02+  =02+0+  *,,2.92E18»-1.11,2.8E4,,7.8E11,0.5/ 

•N2  +  N+  •  N2+  ♦  N  M,1,0,2.02E11,0.81.1.3E4,.7.8E11,0.5/ 

•N  +N  «  N2+  +  E-  •,1,1.0.1.4E13.0.0.6.78E4,,1.5E22,-1.5 
•02  +  M2  ■  NO  +  NQ+  +  E-  ' ,1,1.0.1-38E20,-1.84,1.41E5, ,1 .0E24,-2.S/ 
•NO  +  02*  ND+  +  E-  +  02  ’,1,1 ,0,2.2E1S,-0,3S,1,08E5, ,2.2E26,-2.6/ 
•NO  +  N2-  ND+  +  E-  +  N2  ’ ,1,1 ,0,2.2E15,-0.35,1.08E5, ,2.2E26,-2.5/ 

•0  +  ND+  =  NO  +  D+  ’ ,1,1, 0,3. 63E16, -0.6,5. 08E4. ,1.5E13/ 

•N2  +  0+  =  0  +  B2+  ’,3.4E19,-2.0.-2.3E4,,2.48E19,-2.2/ 

•N  +  N0+  =  NO  ♦  K+  ’,1.0E19,-0.93,6.1E4,.4.8E14/ 

•02  +  M0+  =  NO  +  02+  •,1.8E15,0.17.3.3E4.,1.8E13,0.S/ 

•0  +  N0+=  02  +  N+  ’,1, 1,0, l.34E13,0. 31,7, 727E4,,l.0E14/ 

•END  ’/ 
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F-2.  HYDROGEN  AIR  COMBUSTION 


‘ c . cm/mole , jachimowski ’ /(Ref.  5 2) 

’H2  +  02  =  OH  +  OH  *,l,l,0.1.7E13,,24170./ 

'OH  4  H2  =  H20  +  H  M , 1 ,0.2. 2E13. ,2593./ 

’H  +  02  -  OH  +  H  M,l,0,2.2E14,,e459./ 

»D  +  H2  -  OH  +  H  M,1,0,1.8E10, 1.0, 4481./ 

»0H  +  OH  =  H2D  +  0  M.1.0,6.3E12,,  S49./ 

»H  +  OH  +  H  =  H20  +  H  '.1,1,0,2.2822,-2.0/ 
'H20',6.0/ 

•H+D+M=OH+M  '.l,l,0,6.0E16,-0.6/ 
•H20',S.0/ 

'H+H+H=H2+M  ',1,1,0,6.40817,-1.0/ 
'H2',2.0.'H20’.6.0/ 

'H  ♦  02  4  M  -  HQ2  +  M  ' ,1,1,0,1.7E16. .-503./ 
'H2',2.0,»H20',16.0/ 

'H02  4  H  =  H2  4  02  ' ,1 ,1,0,1, 3E13./ 

*H02  4  H  =  OH  4  OH  ',1,1,0,1.4E14.,544./ 

'H02  4  0  =  OH  4  02  ',1,1.0,1.50813,0.0,478./ 

'H02  4  OH  -  H20  4  02  ',1.1,0,8.0812./ 

'H02  4  K02  =  H202  4  02  ', 1,1,0,2.0812/ 

'H  4  H202  =  H2  4  H02',1,1,0,1.4E12,,1813./ 

'0  4  H202  =  OH  4  H02  M,1,0,1,4E13,  ,3222./ 

'OH  4  H202  »  H2D  4  HD2  '  ,1, 1,0.6. 1E12, ,720./ 

'H  4  H202  =  OH  4  OH  4  M  ' ,1,1, 0,1. 2E17, ,22910./ 
'H2G’.16.0/ 

'0404M=024M  ’,1,1,0,6.0813, ,-503./ 

/ 

'H  4N  4  H  «  N2  4  M  ’,l,l,0,2.8E17.-0.75/ 

/ 

'H  4  02  =  HO  4  0  », 1,1,0, 6. 4E9, 1.0,3172./ 

'M  4  irO  »  M2  4  0  ',1,1,0,1.6813/ 

'M  4  OH  =  NO  4  H  ',1.1,0,6.3E11,0.6/ 

'H  4  HO  4  M  =  HHO  4M  ' ,1,1,0.5.4815, .-302./ 

/ 

'H  4  HHO  -  HO  4  H2  ' ,1,1,0,4.8E12/ 

'0  4  HNO  =  HO  4H20  ’  ,1 ,1,0,5. OEll, 0.5/ 

'OH  4  HMD  =  NO  4  H20  ',1,1,0,3.6813/ 

’H02  4  HHO  =  NO  4  H202  ',1,1,0,2.0512/ 

'H02  4  HO  =  M02  40H  ' ,1,1,0,3.43E12, ,-131./ 

'H  4  H02  -  HO  4  OH  ’,1,1, 0,3. 5814,, 766./ 

'0  4  N02  =  NO  4  02  ',1,1,0,1.0513.. 302./ 

'M  4  ND2  -  NO  4  0  4  H  ', 1 , 1 ,0,1 . 16816, ,33230./ 

/ 

'END'/ 
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Apart  from  the  first  record,  which  is  used  to  identify  its  source  and  date,  the  material 
properties  database  is  divided  into  four  blocks  which  are  read  by  Subroutines  NPCHEMIN, 
NPDIFFIN,  NPVISCIN,  and  NPLAXIN,  respectively.  Table  G-1  is  an  excerpt  from  the 
thermochemical  section  showing  the  data  for  N2.  The  first  line  of  this  block,  which  is  read 
with  the  format 

Table  G-1.  Thermochemical  Data  for  N2 

N2  N  2  .  0  .  0  .  0  .  28.0134  2  4 

3.35324E'^03  .97590E401  .15510E+02  .75  E-19  .S5000B-23  -1.0  E-28 

3.62100E-10  4.0  -15.40393  -.14374E-06  .7S216E-1S- . 14186E-23  .90669E-33 

6  300.0  1000.0  6000.0  15000.0  2S000.0  30000.0 

0.34999E-f01  .10160E-04  .30894E+01 

0.36748E+01-.12081E-02  .23240E-05-.63218E-09-.22577E-12  -17.17441  .23683E+01 

0.28963E+01  . 15155E-02-.57235E-06  .99807E-10-.65224E-14  138.1293  .61615E+01 

0.37270E-^01  .46840E-03-.11400E-06  .11S40E-10-.32930E-15-  32.47182  .13128E+01 

0.96377E+01-.25728E-02  .33020E-06-.1431SE-10  .20333E-15  51.19068  -.37586E+02 
-.51681E+01  .23337E-02-.12953E-06  .27872E-11-.21360E-16  62.96161  ,66217E+02 

1.0  3.0  6.0  6.0 

0.0  .5925684-09  .70954E-t-09  .71042E-409 


(Ix,a8,2x,4(a3,f3.0),f8.4,3i5),  starts  with  the  name  of  the  substance  and  its  chemical  formula. 
Provision  is  made  for  four  element  names.  Also  on  the  first  line  are  the  species  molecular 
weight  {kg/kg-moie),  the  number  of  rotational  degrees  of  freedom,  and  the  number  of 
electronic  energy  levels  (NLEV).  On  the  second  line,  which  is  read  with  format  (lx,7ell.5), 
are  values  for  the  characteristic  vibrational  energy  0  {K),  for  the  standard  heat  of  formation 
(J/kg-mofe)  at  0.0  K  and  1  atm,  for  the  dissociation  energy,  and  for  the  ionization  potential 
(both  in  eV).  The  remaining  three  numbers  on  this  line  are  constants  in  a  polynomial  fit  to 
the  electron-neutral  energy  exchange  cross-section.  On  the  third  line  are  values  for  the  Lenard- 
Jones  collision  diameter,  the  rotational  colbsion  number,  and  Fit  coefficients  for  the  electronic- 
vibrational  energy  relaxation  time.  The  number  of  thermo-fit  segments  (NRNGE)  and  the 
junction  temperatures  are  on  the  next  line.  The  values  of  the  seven  thermo-fit  coefficients 
for  each  range  fill  the  next  NRNGE  lines.  The  last  two  lines  of  data  for  a  species  are  the 
electronic  level  degeneracies  and  corresponding  energies  in  J/kg-mole,  respectively. 

The  boundary  between  the  thermodynamic  data  and  the  diffusion  data  is  marked  by  a 
line  with  the  single  entry  END. 
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An  excerpt  from  the  diffusion  database  which  is  read  by  Subroutine  NPDIFFIN  is  shown 
in  Table  G-2.  This  block  gives,  for  each  species,  the  constants  necessary  to  calculate,  from 
Eq.  (66),  the  binary  diffusion  coefficient  with  the  other  named  species.  These  constants  are 
listed  in  inverse  order,  that  is,  the  coefficient  of  the  highest  power  of  /n  T  is  listed  first  and 
the  constant  term  last.  In  order  to  assure  completeness,  the  symmetry  condition  "Djj  =  "Dj 
is  enforced.  The  diffusion  data  for  different  species  are  separated  by  a  line  with  the  single 
entry  LAST  while  the  boundary  between  the  diffusion  data  and  the  viscosity  data  is  marked 
by  a  line  with  the  single  entry  END. 

The  first  few  lines  of  the  viscosity  data,  which  are  read  by  NPVISCIN,  are  shown  in  Table 
G-3.  The  viscosity  block  begins  with  a  single  line  containing  the  number  of  species  for  which 
data  are  available.  The  following  data  are  read  with  the  format  from  the  second  line.  There 
are  two  lines  of  data  for  each  species:  on  the  first  line,  the  species  name  is  followed  by  an 
integer,  lONIN,  which  connects  a  molecule  or  atom  with  the  ionized  form  of  the  same  species. 
For  instance,  lONIN  =  4  for  NO  and  -  4  for  NO  + .  These  links  are  used  in  setting  the 
species  parameters  for  use  by  NPSUTTON.  As  with  the  diffusion  data,  the  fit  coefficients 
on  the  second  line  are  listed  in  inverse  order  from  their  use  in  Eq.  (65).  The  last  two  entries 
in  each  line  are  the  constants  in  the  Sutherland  viscosity  expression  and  are  no  longer  used. 

The  boundary  between  the  viscosity  data  and  the  relaxation  data  is  marked  by  a  line  with 
the  single  entry  END. 

The  final  block  of  data  in  the  material  properties  database,  as  described  in  Sec.  3.4.1, 
is  read  by  Subroutine  NPLAXIN  and  consists  of  the  parameters  taua,...,taud  in  Eq.  (43) 
of  Sec.  2.5.1,  namely 

p>^tau(j,k}  «  taua*T*naub*«zpCtauc/T**taud) 

for  the  relaxation  time  of  species  j  in  a  bath  of  species  k.  As  shown  in  Table  G-4  ,  there 
is  one  line  for  each  pair  of  species. 
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Table  G-2.  Binary  Diffusion  Data  for  CO 


H2 

-0.491000E~02 

0.143632E+00  -0 

02 

-0.280961E-02 

0.889657E-01  -0 

N2 

<0.292013E-02 

0.915470E-01  -0 

H 

-0.412427E-02 

0.133344E400  -0 

0 

-0.33450&E-02 

0.102934E400  -0 

N 

-0.300961E-02 

0.926432E-01  -0 

OH 

-0.329578E-02 

0.101480E400  -0 

B02 

-0.281176E-02 

0.890325E>01  -0 

K20 

0.368373E-03 

0.728311E-02  -0 

ND 

-0.292217E-02 

0.916079E-01  -0 

AR 

~0.301056E-02 

0.960765E-01  -0 

CO 

-0.293217E-02 

0.919203E-01  -0 

C02 

-0.237574E-02 

0.813037E-01  -0 

cir 

-0.279938E-02 

0,866140E-‘01  -0 

HCN 

-0.974226E-03 

0.448550E-01  -0 

N20 

-0.243124E-02 

0.825726E-01  -0 

N02 

-0.287876E-02 

0.945216E-01  >0 

C 

-0.310294E-'02 

0.964079E-0i  -0 

C2 

LAST 

-0.293334E-02 

0.919424E-01  -0 

.  lS7448E-t'01  0 . 932089E+01  -0 . 236352E-K}2 
.  105720E+01  0 . 724306E+01  -0 . 220271E+02 
.  107718E+01  0 . 729354E+01  -0 . 220264E+02 
.161427E+01  0.103288E+02  -0.265753E+02 
.118857E+01  0.77593EE+01  -0 . 222797E+02 
.107007E+01  0.71S20BE+01  -0.212231E+02 
.117263E+01  0.768092E+01  -0.221662E+02 
.105797E+01  0.724703E+01  -0.220419E+02 
.32363iE-i>00  0.46e910E+01  -0 . 194570E+02 
.107787E+01  0.729693E+01  -0.220496E+02 
.115U9E+01  0.779324E+01  -0.232613E+02 
.1081S6E+01  0,731631E+01  -0.220796E+02 
.104026E+01  0.755574E+01  -0.236811E+02 
.100572E4D1  0.684965E401  -0.210421E+02 
.708852E400  0.638614E401  -0 . 22&367E402 
.  104943E401  0 . 757232E401  -0 . 236628E+02 
.116502E401  0.804363E401  -0 . 24227&E402 
.110070E+01  0.730243E+01  -0.2144S3E+02 
.108161E+01  0.731B62E+01  -0.220272E+02 
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Table  G-3.  Viscosity  Data 


41 

(7E10.4) 

02  1 


0.0000 

N2 

2 

0.00 

0.0484 

-0.1455 

-8.9231 

1.4584E-6 

110.33 

0.00 

M 

5 

0.0 

0.0203 

.4329 

-11.8153 

1.4S84E-6 

110.33 

0.00 

Q 

3 

0.0 

0.0120 

0.5930 

-12.3805 

1.458^-6 

110.33 

0.00 

NO 

4 

0.00 

0.0205 

0.4257 

-11.5803 

1.5864E-6 

110.33 

0.00 

N0+ 

-4 

0.00000 

0.0452 

-0.0609 

-9.4596 

1.4564E-6 

110.33 

0.0913 

-3.3178 

45.1426 

-270.3958 

586.3300 

1.4584E-6 

110.33 

Table  G*4.  Vibrational  Relaxation  Rate  Coefficients 

’M2','M2M.6909E-11.0.0,220.0,0.333 
‘CO',  'ARM. 6986E-11, 0.0, 213., 0.333 
'CO' , 'CO' ,6.2397E-11,0.0.175. .0.333 
'CO* , 'HE' ,1.3127E-09, 0.0.99. 0,0.333 
'CO' , 'H2' ,3.0959E-09.0.0,68.0,0.333 
'02' , 'AR' ,6.0211E-11.0.0,166. .0.333 
'02' , '02' .2.0872E-10,0.0,129. ,0.333 
'02' , 'HE' , 2. 5162E-09, 0.0.67. 0.0. 333 
'02' . 'H2» ,4.78S3E-09,0.0,42.0.0.333 
'F2',’F3',1.3069E-09,0.0.65.0.0.333 
'CL2' , 'CL2' , 1 . 1968E-09. 0.0.58. ,0.333 
’BR2' , »BR2» , 1 . 1625E-09.0.0,48. ,0.333 
'I2*,'I2',2.3233E-09,0.0.127.0,0.333 
'MO', 'jrD’,7.0219E-ll,0.0,168.0,0.333 
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The  material  properties  database  currently  contains  data  for  the  species  listed  in  Table  H-1 . 


Table  H-I.  Tabulated  Species 


E- 

C 

H 

N 

0 

AR 

C2 

H2 

N+ 

0+ 

AR+ 

C2- 

H20 

N2 

0- 

C3 

H20+ 

N2+ 

02 

CH 

H202 

N20 

02+ 

CH2 

H3+ 

NCO 

02- 

CH20 

HCN 

NH 

OH 

CH20H 

HCO 

NH3 

0H+ 

CHS 

HCO+ 

NO 

OH- 

CH30H 

H£ 

NO+ 

CH4 

HE+ 

N02 

CN 

HNCO 

CO 

HNO 

C0+ 

H02 

C02 
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APPENDIX  I 

A  TYPICAL  APPLICATION 

A  typical  application  of  NEQPAK  to  a  hypothetical  problem  is  outlined  below.  There 
are,  of  course,  many  flow  problems,  both  internal  and  external,  for  which  an  understanding 
of  thermochemical  effects  is  important.  There  are  also  many  flow  solvers  now  available.  In 
order  to  be  concrete,  a  particular  class  of  problem  and  a  particular  computational  approach 
will  be  assumed.  Only  the  procedures  which  involve  the  chemistry  or  thermophysical  properties 
will  be  discussed. 


The  problem  to  be  solved  is  that  of  determining  the  electromagnetic  properties  of  the 
plasma  about  a  reentry  body.  An  iterated  implicit  flow  code  with  shock  fitting  is  used.  The 
reaction  set  “blottner”  from  Appendix  F-1  and  thermal  equilibrium  are  selected. 

•  PRELIMINARIES 

—Choose  values  for  INPR,  INREL,  and  INPTH  which  will  not  interfere  with 
other  input-output  units  needed  by  the  code. 

—Open  INPR,  INREL,  and  INPTH  to  read  formatted  files. 

—Set  maximum  array  dimensions  by  editing  file  NPDIMEN  (Table  4). 

•  INPUT 


—Set  NSFS  =2;  SPECFS(1)  =  02.  SPECFS(2)  =  N2, 

CON(l)  =  .22,CON(2)  =  .78 

—Set  MODEL  =  1,  NVIB  =  0,  VIB  =  dummy,  HALT  =. FALSE. 
—CALL  NPCHEMIN( . )  (Table  5). 


— After  return  from  NPCHEMIN: 


•SPECFS(l)  =  E-,  CON(1)  =  1.0E-15 
*SPECFS(2)  =  O,  CON(2)=  1.0E-I5 
*SPECFS(3)  =  N,  CON(3)  =1.0E-15 
•SPECFS(4)  =  O2,CON(4)  =  0.22 
•SPECFS(5)  =  N2,  CON(5)=0.78 
•SPECFS(6)  =  NO-h,  CON(6)=1.0E-15 
•SPECFS(7)  =  NO,  CON(7)=1.0E-1S 
*NSPEC  =  7 


•NTMP  =  1 
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—Set  GAM(IS.1)  =  C0N(IS).  IS  =  1,NSPEC 

—CALL  NPWEIGHT  (GAM.l.l,l,’mole/mole’.WT)  to  get  free-stream 
molecular  weight  WT(1). 

•  POST-SHOCK  CONDITIONS 

—Set  up  initial  shock  geometry  and  free-stream  conditions  dependent  on  this 
geometry.  In  this  hypothetical  case,  a  uniform  free  stream  is  assumed. 

— Since  the  shock  has  infinitesimal  thickness,  there  is  no  time  for  reactions  to 
occur  in  a  gas  parcel  crossing  the  shock.  Therefore,  set  GAM(IS,IP) 
=  CON(IS)forIS  =  I,NSPECandIP  =  1  ,NSH  where  NSH  is  the  number 
of  points  immediately  behind  the  shock. 

-CALL  NPTHERM  (TFS,TFS,TFS,l,NSH,l,MODEL,HSF,HSF,HSF,CPS, 
CPS, CPS)  where  TSF  is  free-stream  temperature  and  HSF  is  the  correspond¬ 
ing  species  enthalpy. 

— Calculate  EGAM(IS,IP)*HSF(IS,IP)  at  each  point  on  the  shock  to  get  the 
static  enthalpy  of  the  free-stream  mixture. 

— From  momentum  and  mass  conservation,  estimate  a  velocity  and  then  a  post¬ 
shock  static  enthalpy  h.  Also  estimate  a  corresponding  temperature  T. 

—CALL  NPHFINDT{/i,GAM,  1 ,1  ,NSH,1  ,T)  to  get  the  post-shock  temperature. 

— Change,  say,  the  post-shock  density  and  repeat  the  last  two  steps  until 
convergence  is  obtained. 

•  INITIALIZE  SHOCK  LAYER 

— Estimate  temperature  T,  pressure  P,  and  mole-fraction  distributions  consistent 
with  shock  and  surface  conditions. 

—CALL  NPWEIGHT  to  get  the  molecular  weight  and  then  calculate  the  density 
RO  from  the  equation  of  state. 

— Save  the  mole  fractions  in  CHI  and  CALL  NPTOCON(GAM,NP,’mole/- 
mole’,RO,WT)  to  convert  GAM  to  concentration  units. 
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— CALL  NPTHERM  to  gel  the  enthalpy  and  specific  heat  distributions, 

—CALL  NPTRANS{LOW,NP,INC,T,T,CPS,CPS,CPS,GAM,D,VIS,EMU, 
TKT,TKT,TKT)  to  calculate  the  transport  properties. 

— Since  an  implicit  algorithm  will  be  employed  by  the  flow  solver,  set 
DERIV  =  .TRUE. 

— Since  only  air  is  involved,  there  will  be  virtually  no  low-temperature  chemical 
activity,  so  to  save  computing  time,  set  TMIN = 800.0. 

—CALL  NPSOURCE(GAM,T.DERIV,LOW,NP,INC,GAIN.LOSS,WDOT, 
DWDTM,DWDG). 

— The  user  now  has  sufficient  information  to  assign  initial  values  to  the  variables 
at  every  point  of  the  shock  layer  and  also  to  initialize  the  integration  algorithm. 

•  SOLVE  CONSERVATION  EQUATIONS 

— After  the  first  pass,  new  values  for  global  density,  species  concentrations, 
and  mixture  static  enthalpy  will  be  available  at  each  point  of  the  shock  layer. 

— Call  NPHFINDT  to  get  the  temperature  at  each  point, 

— Call  NPTHERM  to  get  species  specific  heats. 

— Call  NPTRANS  to  reevaluate  the  transport  properties. 

— Check  magnitude  of  temperature  and  composition  changes.  If  the  changes 
are  small  enough,  set  DERIV  =  .FALSE,  to  save  CPU  time. 

— Call  NPSOURCE  to  get  new  source  terms. 

— Reapply  solver  until  global  convergence  is  obtained. 

•  REPORT  RESULTS 

— Organize  output  data. 

—Call  NPHFINDT  to  get  final  values  for  temperature  T. 
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NOMENCLATURE 

Constant  in  forward  rate  coefficient,  Eq.  (28),  m^/kg — mole  or 
m^/kg — mol^ 

Aij  Coefficient  defined  by  Eq,  (47), 

Or  Effective  activation  temperature  aC„  Sec.  2.4.2,  K 

Ojs  Thermo-fit  coefficient,  Eqs.  (62)  —  (64), 

Constant  in  forward  rate  coefficient,  Eq.  (28) 

Coefficient  in  viscosity  fit,  Eq.  (65) 

Aaivation  temperature,  Eq.  (28),  K 

Cp  s  Constant  pressure  molar  heat  capacity  of  species  j,  J/kg  —  mole/K 

Df  Constant  in  reverse  rate  coefficient,  Eq.  (67) 

Dg  Ambipolar  diffusion  coefficient,  mVj 

Dg  Effective  diffusion  coefficient,  Eq.  (16),  m^/s 

DJ  Thermal  diffusion  constant,  /«"* 

Binary  diffusion  constant,  m^/s 
Diffusion  vector  defined  in  Eq.  (14) 
djf  jg  Coefficient  in  Eq.  (66) 

E  Total  specific  energy  of  mixture,  J/kg 

Ef  Constant  in  reverse  rate  coefficient,  Eq.  (67) 

e  Magnitude  of  charge  on  electron,  C 

^  Mixture  internal  energy  of  type  /,  J/kg 
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4  Species  internal  energy  of  type,  /,  J/kg 

Constant  in  reverse  rate  coefficient,  Eq.  (67) 

AF  Change  in  free  energy,  J/kg  —  mole 

ith  component  of  body  force  on  species  s,  N/kg  —  mole 
Gibbs  free  energy  of  species  s,  J/kg  —  mole 
Qs  Rate  of  appearance  of  species  s,  kg  —  mole/m^/s 

Or.s  Contribution  of  reaction  r  to  kg  —  mole/m^/s 

Q]  Degeneracy  of  electronic  state  / 

Pseudo-temperature,  Eq.  (39) 

H  Total  enthalpy  of  mixture,  J/kg 

A 

hg  Mass-specific  static  enthalpy  of  species  s,  J/kg 

hg  Static  enthalpy  of  species  s,  J/kg  —  mole 

Ao,f  Heat  of  formation  at  0  K  and  1  atm  pressure  of  species  s,  J/kg  —  mole 

Translation-rotation  enthalpy  of  species  s,  J/kg  —  mole 
hg  Vibration-electronic  enthalpy  of  species  s,  J/kg  —  mole 

h*  Vibrational  enthalpy  of  species  s,  J/kg  —  mole 

hi  Electronic  enthalpy  of  species  J,  J/kg  —  mole 

Is  Ionization  energy  of  species  s.  Sec.  2.6,  J/kg  —  mole 

Jfj  Element  of  Jacobian 

Kl  Equilibrium  constant,  Eq.  (30) 
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k/f  Forward  rate  coefficient  of  reaction  r,  n^/kg  —  mole  or  m^/kg  —  mol^ 
kj.  Reverse  rate  coefficient  of  reaction  r,  m^/kg  —  mole  or  m^/kg  —  mol^ 

L  Debye  length,  m 

Cg  Rate  of  disappearance  of  species  s,  kg  —  mole/m^/s 

Contribution  of  reaction  r  toC^,  kg  —  mole/m^/s 
M  A  generic  molecule,  kg  —  mole 

M's  Molecular  weight  of  species,  kg/kg  —  mole 

mg  Electron  mass,  kg 

Avogadro's  number 

N  LEV  Number  of  electronic  energy  levels 

Ng  Number  of  reactions 

iVj  Number  of  species 

fig  Electron  number  density,  m~^ 

Nj-  Number  of  temperature  types  ^ 

p  Static  pressure  of  mixture,  N/m^ 

Partial  pressure  of  species  s,  N/m^ 

Vij  Probability  of  vibrational  energy  exchange 

Q  Radiant  energy  transfer  rate,  J/m^/s 

g  Heat  flux  vector,  J/nP-fs 

H  Universal  gas  constant,  J/kg  —  mole/K 
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5  Electrostatic  shielding  ratio 

Electron  energy  loss  rate  by  impact  ionization,  J/m^/s 

Electron  energy  loss  rate  by  collision,  J/m^/s 

Vibration  energy  gain  by  chemical  reaction,  J/m^/s 

Vibration  energy  gain  by  exchange,  J/rn^/s 

Vibrational  energy  gained  by  electron  impact  for  species  s,  W/rr^ 

Vibrational  energy  source  for  species  s,  W/m^ 

Dissociation  temperature  of  species  s,  K 
Generic  temperature,  K 

Tg  Translational  temperature  of  free  electrons,  K 

Tf  Generic  internal  temperature,  K 

T  Translational  temperature,  K 

Ty  Temperature  of  all  internal  modes  other  than  rotation,  K 

Ty  Vibrational  temperature,  K 

T  Pseudo-temperature,  Eq.  (41),  K 

Tq  j.  Pseudo-temperature,  Eq.  (40),  K 

t  Time,  s 

T  Natural  logarithm  of  T 

U  Mass  averaged  flow  velocity,  m/s 
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MS 


Z 


Z 

a 


0 

r 

rteJr) 


A, 

«r 


e/ 


Ah  component  of  U,  m/s 

Diffusion  velocity  of  species  s,  m/s 

Ah  component  of  m/s 

Vibrational  quantum  number 

Rate  of  production  of  species  s,  kg  —  mole/m^/s 

Mass  fraction  of  species  s 

Vibrational  partition  function,  £q.  (35) 

Electronic  partition  function 

Coefficient  defined  by  Ref.  10 
Number  of  atoms  of  element  i  in  species  j 
Exponent  of  Ty  'm  Park’s  average  T^,  Sec.  2.4 
Factor  in  Eq.  (34) 

Total  concentration,  kg  —  mole/rn^ 

Incomplete  gamma  function 
Species  concentration,  kg  —  mofe/rn^ 

Dissociation  energy  of  species  s,  J/kg  —  mote 
Change  in  number  of  moles  across  reaction  r 
Kronecker  delta;  =  1  if  /  =  j,  =0  if  /  ^  j 
Energy  of  Ah  electronic  state,  J/kg  —  mole 
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ij  Viscosity  of  mixture,  Ns/rri^ 

Characteristic  vibrational  temperature  of  species  s,  K 
X  Boltzmann  constant,  J/ particle/ K 

\  Frozen  thermal  conductivity  of  species  s,  J/m/s/K 

Type-/  thermal  conductivity  of  mixture,  J/m/s/K 
A  Electrostatic  shielding  function 

Stoichiometric  coefficient  of  reactants  in  reaction  r 
p’rj  Stoichiometric  coefHcient  of  products  of  reaction  r 

Ps  Characteristic  vibrational  frequency, 

Q  Density  of  mixture,  kg/m^ 

Gj  Density  of  species  s,  kg/m^ 

Tj  Vibrational  relaxation  time,  Eq.  (44),  s 

Probability  of  dissociation,  Eqs.  (39)  -  (41) 

Xg  Mole  fraction  of  species  s 

'Jy  Forward  coupling  factor,  Sec.  2.4.2 

Reverse  coupling  factor.  Sec.  2.4.2 
l/’;  Dimensionless  energy  of  /th  electronic  level 

Rate  of  Increase  of  internal  energy  of  type  /,  J/kg/s 
u-ij  Reduced  mass  of  collision  partners  /  and  j,  kg 

^5  Mobility  of  ion  5,  m^/s/voU 
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